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Abstract 


This  report  summarizes  the  main  results  of  research  conducted  during  the  summer  of  2007  on 
tasking  a  finite  number  of  cooperative  agents  to  randomly  emerging  targets  for  their  removal. 
Faults  occur  when  some  agents  engaged  in  a  mission  are  expired.  Agents  are  subject  to  threat  at  a 
level  determined  by  the  number  of  targets  present.  On  the  other  hand,  the  rate  at  which  a  target  is 
removed  depends  on  the  number  of  cooperative  agents  assigned  to  it.  Faults  effectively  change  the 
network  architecture  and  therefore  degrade  the  network  performance.  Designs  of  control  policies 
that  determine  the  number  of  agents  assigned  are  based  on  the  network  life  when  expired  agents 
cannot  be  replenished,  and  on  the  network  availability  when  expired  agents  are  replenished  at  a 
certain  rate.  Tasking  process  is  described  by  a  discrete  event  system  in  the  form  of  a  queuing 
network,  where  agents  are  servers  and  targets  are  customers.  Optimal  policies  are  determined  by 
solving  a  Markov  decision  problem.  To  facilitate  the  readers  understanding  of  the  motivation,  and 
of  the  problem,  the  agents  are  specialized  to  networked  pairs  of  airborne  sensors  that  are  tasked  to 
locate  non-cooperating  microwave  transmitters  as  targets. 

This  work  has  resulted  in  a  Master’s  thesis  by  Yan  Guo  whose  graduate  study  was  partially  sup¬ 
ported  by  this  grant,  a  paper  submitted  to  2008  IFAC  World  Congress,  and  another  paper  submitted 
to  International  Journal  of  Control,  Automation,  and  Systems. 

Some  simulation  results  are  also  included.  These  results  are  not  fully  satisfactory  and  are 
incomplete.  Continued  work  is  desired. 


1  Introduction 


This  work  considers  tasking  a  finite  number  of  cooperative  agents  to  randomly  emerging  targets 
for  their  removal.  It  focuses  on  enhancing  the  network  performance  in  the  face  of  expiration  of 
its  agents.  The  resulting  network  is  said  to  be  fault-tolerant.  When  the  agents  considered  are  net¬ 
worked  airborne  sensors,  the  mobility  and  a  multiplicity  of  the  sensing  nodes  make  fault-tolerance 
possible.  Fault-tolerant  tasking  in  this  work  is  achieved  by  implementing  operation  policies  opti¬ 
mized  for  network  availability. 

Control  of  networked  multiple  agents  has  been  an  intensively  discussed  topic  recently  in  the 
controls  literature  [1],  [12],  for  example,  describes  a  pursuit  evasion  game,  where  mobile  agents 
are  to  chase  and  capture  multiple  moving  targets  in  a  minimum  amount  of  time,  and  a  network  of 
stationary  sensors  serves  to  help  enhance  the  target  observability  in  the  game. 

With  unmanned  aerial  vehicles  (UAVs)  replacing  stationary  networks  and  manned  vehicles, 
significant  improvements  in  network  performance  can  be  expected.  Networking  in  a  hostile  envi¬ 
ronment,  however,  poses  new  challenges.  Data  exchange  inherent  to  a  networked  operation  and 
prolonged  mission  time  due  to  poor  execution  expose  the  otherwise  passive  location  sensors,  thus 
increase  the  likelihood  of  the  vehicles  being  destroyed. 

Examine  a  situation  where  the  motion  of  two  unmanned  aerial  vehicles  (UAV)  and  a  hostile 
radar  lie  within  a  plane,  as  illustrated  in  Fig.  1.1.  Let  us  assume  that  the  two  vehicles  are  equipped 
to  acquire  both  the  time  difference  of  arrival  and  the  frequency  difference  of  arrival  of  the  radar 
signal  [6].  The  sensors  are  mounted  on  the  vehicles,  and  thus  are  subject  to  the  same  speed  and 
curvature  constraints  as  that  of  the  vehicles.  The  sensors  are  passive  nodes  when  acquiring  data 
from  the  transmitter,  but  become  active  when  exchanging  data  between  them  in  order  to  provide  a 
location  estimate. 

It  can  be  seen  from  Figure  1  that  at  least  two  sensor  carrying  vehicles,  which  make  both  a 
TDOA  measurement  and  an  FDOA  measurement,  are  needed  in  a  2-dimensional  setting  to  locate 
the  target.  A  noiseless  measurement  of  time  difference  of  arrival  by  a  pair  of  sensors  on  the  two 
vehicles  is  given  by 


St  =  ~W{x2  —  xe)2  +  (y2  ~ye )2 

-  V(xi  ~xe)2  +  (2/i  —  2/e)2],  (1.1) 

and  a  noiseless  measurement  of  frequency  difference  of  arrival  by  the  same  pair  of  sensors  is  given 
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Figure  1.1:  Transmitter  location  using  a  TDOA  measurement  and  an  FDOA  measurement  by  two 
airborne  sensors. 


by 


feA'X 2  -  Xe)u2  +  (?/2  ~  Ve)V2 


c  \J (x2  ~  xe)2  +  (y2  -  ye )2 


(1.2) 


y/(xi  —  Xe)2  +  (yi  —  |/e)2 


where  (xe,ye)  is  the  transmitter  location  to  be  estimated,  (x\  .  y\  )  and  (x2,y2)  are  the  positions 
of  the  two  vehicles,  respectively,  (ui,Vi)  and  (u2,v2)  are  the  velocities  of  the  vehicles,  fe  is  the 
carrier  frequency  of  the  transmitted  signal,  and  c  is  the  speed  of  light. 

Since  the  measurements  are  always  noisy,  multiple  measurements  are  needed  for  an  accurate 
location  estimation  of  the  emitter.  Such  measurements  can  be  distributed  temporally  along  the 
trajectories  of  motion  of  a  pair  of  sensors,  or  spatially  over  multiple  pairs  of  sensors,  or  both. 
Measurements  made  by  multiple  pairs  of  sensors,  which  form  a  network,  offer  greater  degree  of 
fault-tolerance,  and  greater  potential  for  improved  speed  and  accuracy  in  target  location.  The 
reader  is  referred  to  [6]  and  [14]  for  more  detailed  discussion  on  methods  for  location  estimation 
and  accuracy  analysis. 

A  tasking  problem  that  is  specific  to  this  application  refers  to  that  of  allocating  a  finite  number 
of  sensor  pairs  to  randomly  emerging  microwave  transmitters  to  maximize  the  network  availability. 
A  tasking  policy  that  is  too  greedy  tends  to  exhaust  resources  before  the  arrival  of  unanticipated 
radars,  whereas  a  tasking  policy  that  is  too  conservative  tends  to  lengthen  the  exposure  of  the 
sensor  carrying  vehicles.  Tasking  is  treated  as  a  server  allocation  problem  of  a  queuing  network. 
Optimal  policies  are  sought  as  the  solutions  of  Markov  decision  problems. 

The  report  is  organized  as  follows.  Chapter  2  modeling  the  tasking  process  for  a  small  scale 
sensor  network.  Chapter  3  designs  supervisory  control  policies  for  optimal  tasking  for  the  cases 
where  lost  sensors  can  and  cannot  be  replenished  by  solving  appropriate  Markov  decision  prob¬ 
lems.  Chapter  4  evaluates  the  network  performance  in  terms  of  expected  network  life  and  steady- 
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state  availability.  Chapter  5  concludes  the  report. 
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2  Modeling  of  tasking  process 


An  optimized  tasking  is  one  that  maximizes  the  expected  life  of  the  network  where  the  lost  air¬ 
borne  sensors  cannot  be  replenished,  or  one  that  maximizes  the  expected  steady-state  availability 
of  the  network  where  the  lost  sensors  can  be  replenished.  In  this  study  fault-tolerance  refers  to  the 
network’s  tolerance  to  vehicle  loss. 

Figure  2.1  is  a  queuing  network  model  of  a  six-sensor,  finite  target  population  tasking  process, 
where  each  server  represents  a  pair  of  sensors  capable  of  independently  locating  a  target  to  a  certain 
accuracy  in  the  absence  other  pairs,  and  a  customer  is  a  randomly  emerging  target. 

Each  customer  resides  in  the  queue  or  a  server  is  regarded  as  a  detected  target  which  is  being 
or  to  be  served  by  one  or  more  servers  or  sensor-pairs.  Service  is  complete  as  soon  as  the  target 
location  is  determined  to  a  required  accuracy.  A  target  is  then  considered  removed.  A  sensor-pair 
allocated  to  a  target  is  tied  to  the  target  until  its  service  is  complete,  or  the  life  of  the  sensor-pair 
is  terminated,  whichever  comes  first.  The  three  delay  elements  of  average  delay  1/A  imply  that 
target  population  is  limited  by  three  at  any  given  time.  A  new  target  is  generated  or  replenished  at 
a  delay  element  with  rate  A  upon  the  service  completion  of  a  target  at  one  or  multiple  servers. 

An  supervisory  control  policy  determines  whether  to  allocate  one,  or  two,  or  three  pairs  of 
sensors  to  each  reported  target,  with  a  corresponding  mean  service  time  of  1/ni,  (>)1/ /z2,  (> 
)l//i3,  respectively,  where  //,  denotes  the  service  rate  of  committing  i  pairs  of  sensors  to  a  target. 
Given  the  sensing  mechanism,  the  mean  service  time  by  a  single  pair  of  sensors  is  in  the  range  of 
seconds  to  tens  of  seconds,  dominated  by  the  time  required  to  adjust  sensor  positions  and  velocities 
for  continued  data  collection,  exchange,  and  processing  needed  for  target  location  to  a  required 
accuracy.  Each  sensor-pair  has  a  mean  lifetime  l/u0  >  \/v\  >  l/u2  >  1  /z/3,  depending  on  the 
threat  level  quantified  by  the  number  of  targets  present  as  indexed  by  the  subscript.  l/z/0  is  the 
server  life  representing  the  expected  natural  endurance  of  a  vehicle,  which  is  “often  an  hour  or  so 
at  best”  [13].  It  also  reflects  sensor  lives  affected  by  undetected  targets.  The  network  is  said  to 
have  expired  when  there  is  no  longer  a  single  surviving  sensor-pair. 

Tasking  process  model  is  built  in  this  study  with  the  premise  that  event  life  distributions  have 
been  established  for  the  process  of  target  arrival  (exp(A)  =  1  —  e~xt),  the  process  of  target  location 
(exp(/ij)),  the  process  of  loss  of  a  sensor-pair  (exp(z/j)),  and  the  process  of  sensor  replenishment 
(exp(cu))  when  new  sensor  carrying  vehicles  are  supplied  for  an  expired  network.  Since  all  event 
lives  are  assumed  to  be  exponentially  distributed,  the  database  unit  can  be  conveniently  modeled 
as  a  Markov  chain  specified  by  a  state  space  X ,  an  initial  state  probability  mass  function  (pmf) 
ttx  (0) ,  and  a  set  of  state  transition  rates  A,  /q,  i and  u. 

A  state  name  is  coded  with  a  4-digit  number  indicative  of  the  number  of  targets  present  and 
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(c)  Three  sensor-pair/target  allocation 

Figure  2. 1 :  A  queuing  network  model  of  a  three  sensor-pair,  finite  target  population  airborne  sensor 
network. 
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the  network  configuration.  A  valid  state  representation  is  given  by  QS,  where  queue  length  0  e 
{0,1,  2,  3},  and  server  state  S'  =  (i,  j,  k),  with  i  e  {0,1,  2, 3, 4},  and  j  e  {0, 1, 2,  3, 4,  5},  and 
k  G  {0, 1,  3, 4, 5}.  A  server  state  “0”,  represented  by  the  value  of  i,  or  j,  or  k,  indicates  an  idle 
sensor-pair,  a  “1”  indicates  a  target’s  being  located  by  one  server  (or  one  sensor-pair),  a  “2”  and 
a  “3”  indicate  that  a  target’s  being  located  by  two  and  three  cooperating  servers,  respectively,  a 
“4”  indicates  a  lost  server,  and  a  “5”  indicates  that  the  lost  server  has  been  tied  to  another  server 
in  serving  a  target.  The  expired  network  requires  4  distinct  states  to  memorize  the  possible  queue 
length  distributions.  Note  that  this  state  specification  has  assumed  homogeneous  sensor-pairs  and 
homogeneous  targets,  and  has  made  use  of  the  symmetry  which  results  in  37  states.  A  set  of 
alternative  state  names  are  assigned  from  X  =  {1, 2, ...,  37}  with  0000  mapped  to  x  =  1  and  the 
network  expiration  states  mapped  to  x  =  34,  35, 36,  and  37. 

Events  that  trigger  the  transitions  and  the  corresponding  transition  rates  are  given  as  follows. 
An  emerging  target  enters  with  rate  (3  —  Q)  x  A.  A  target  is  located  by  one  sensor-pair  with  rate 
/ii,  and  i(>  1)  cooperative  sensor-pairs  with  rate  /p.  In  the  latter  case,  the  i  servers  are  configured 
as  a  single  hyper-exponential  server  with  i  parallel  stages  [4].  An  arriving  target  enters  any  one 
of  the  servers  with  probability  1/i,  which  has  a  service  time  distribution  exp (  //,;  ).  When  service  is 
completed,  the  target  is  removed,  while  no  new  target  can  enter  service  when  the  hyper-exponential 
server  is  busy.  The  service  time  distribution  of  a  hyper-exponential  server  is 


W  =  E7(1- 


=  1  -e 


—Hit 


3= 1 


(2.1) 


which  assumes  homogeneity  of  the  servers.  Loss  of  a  sensor-pair  occurs  at  rate  mu 0  when  the 
network  is  idle  with  m  remaining  sensor-pairs,  mul  when  one  target  emerges,  and  mul  when 
i(>  1)  targets  emerge.  Replenishment  process  begins  at  the  network  expiration  with  rate  u>.  If  one 
of  the  sensor-pairs  is  lost  while  locating  a  target  with  other  sensor-pairs,  the  surviving  sensor-pairs 
continue  to  locate  the  target  at  the  same  rate.  This  is  a  simple  way  to  memorize  the  service  already 
being  provided  without  resorting  to  a  more  complex  model. 

Let  X  e  X  denote  the  random  state  variable  at  time  t.  The  set  of  state  transition  functions  is 
given  by 

Pi, jit)  =  P[X{t)  =j\X(0)  =  i],  i,j  =  1,2,..., 37.  (2.2) 


The  continuous-time  Markov  chain  can  be  solved  from  the  forward  Chapman-Kolmogorov  equa¬ 
tion  [4],  [9] 

Pit)  =  P{t)Q{u(x)),  P{  0)  =  /,  Pp)  =  \pitj{t)\  (2.3) 

and  Q(u(x))  is  called  an  infinitesimal  generator  or  a  rate  transition  matrix  whose  (i.  j)th  entry 
is  given  by  the  rate  associated  with  the  transition  from  current  state  i  to  next  state  j.  Table  1 
summarizes  information  contained  in  transition  rate  matrix  Q(u(x)).  Control  variable  u(x)  will  be 
discussed  shortly.  State  probability  mass  function  at  time  t 


VT it)  =  [7Ti it)  7 r2(t)  ■  ■  ■  TT 37(f)],  t  >  0  (2.4) 

can  be  solved  from 

7i(£)  =  n(t)Q(u(x)),  given  n(t  =  0).  (2.5) 
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A  Markov  chain  for  the  tasking  process  of  Figure  2  has  been  established  so  far.  Since  transition  rate 
matrix  Q  is  dependent  on  control  actions,  the  state  transition  functions  Pij(t )  are  being  controlled, 
and  so  are  the  state  probabilities.  Rate  transition  matrix  Q  is  given  in  the  form  of  a  table  in  Table 
1. 
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Figure  2.2:  Transitions  and  transition  rates  of  the  tasking  process 
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3  Design  of  supervisory  control  policy 


Several  possible  supervisory  control  policies  associated  with  tasking  are  examined.  An  aggressive 
policy  allocates  as  many  available  sensor-pairs  to  as  many  targets  present;  A  greedy  policy  allocates 
all  available  sensor-pairs  to  one  target  at  a  time;  A  conservative  policy  always  allocates  only  one 
sensor-pair  to  every  target  present  to  reserve  assets  in  anticipation  of  new  targets.  In  addition,  four 
optimal  policies  have  been  attempted  to  minimize  the  cost  of  sensor  loss,  threat  level,  unattended 
targets,  and  time  needed  to  replenish  upon  network  expiration,  respectively. 

The  optimal  policies  are  obtained  by  solving  Markov  decision  problems  of  appropriate  penalty 
functions.  A  discrete-time  Markov  chain  model  suitable  for  this  purpose  can  be  derived  under  each 
cost  criterion  by  the  application  of  a  uniformization  procedure  [9] 

7t(4+i)  =  7T  (tk)[I  +  ^Q(u(xk))\,  (3.1) 

where  the  uniform  rate  p  is  greater  than  any  total  outgoing  transition  rates  at  any  states  of  the 
original  continuous-time  Markov  chain  (2.5). 

Each  Markov  decision  problem  considered  in  this  report  assumes  that  a  cost,  denoted  by 
C(i,u),  is  incurred  at  every  state  transition,  where  i  is  the  state  entered  and  u  is  a  control  action 
selected  from  a  set  of  admissible  actions  [4],  [2].  A  solution  amounts  to  determining  a  stationary 
policy  7 r  =  {u(xk),  k  —  0, 1,  •  •  •  }  that  minimizes  the  following  expected  total  discounted  cost 


OO 

Vn(x0)  =  EnJ2  <xkC{Xki  uk)  (3.2) 

k= 0 

where  0  <  a  <  1  is  a  discount  factor.  C(Xk,  uk)  in  each  of  the  four  Markov  decision  problems 
takes  the  form  of  total  number  of  lost  sensor-pairs,  cu-1  at  the  state  of  network  expiration,  u,  at  the 
state  where  it  is  the  server  loss,  and  Q  at  the  state  where  it  is  the  queue  length,  respectively. 

Let  Xk  e  (1,  2,  -  -  -  ,37}  denote  the  random  state  variable  at  tk  =  k/p  in  the  discrete  time 
Markov  chain.  Control  action 

{1,  allocate  one  sensor-pair  to  a  target 

2,  allocate  two  sensor-pairs  to  a  target  (3.3) 

3,  allocate  three  sensor-pairs  to  a  target 

Note  that  the  indicator  functions  in  Table  1  are  defined  follows 


(  1,  u  =  i 

\  0,  otherwise 


A  =  1,2,3. 


(3.4) 
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It  is  known  [4,  2]  that  under  the  condition  0  <  C(j,  u)  <  oo  for  all  j  and  all  u  that  belongs  to 
some  finite  admissible  sets  Uj,  the  minimum  cost  V*(i)  satisfies  the  following  optimality  equation: 


V(i) 


37 


mm 

u&Ui 


C(i ,  u)  +  a  ^  pij V (j)  >  ,  ueUh 


j= i 


(3.5) 


i  —  1,  •  •  •  ,37,  where  pjj  is  the  ( i,j)th  entry  of  I  +  ^(3(u(xfc)). 

The  solution  to  (3.5)  can  be  obtained  via  linear  programming  [3,  2].  In  this  case,  the  set  of 
optimality  equations  is  turned  into  a  set  of  affine  constraints  on  the  set  of  optimization  variables 
{V (/')},  and  the  problem  can  be  formally  stated  as  follows. 


Maximize 

1/(1)  +  1/(2) +  •••  +  1/(36)  +  1/(37) 

(3.6) 

Subject  to 

V(i)  >0,  i  G  X  =  (1,  -  -  -  ,37} 

(3.7) 

V(i)  <  [C{i,u)  +  a^2pitjV(j)}  |u, 

(3.8) 

j 


Vw  eUi,  i  e  X. 

In  the  tasking  process  considered,  Uj  is  nonempty  only  at  state  j  =  1,  2,  4,  7,  13,  14,  15,  19, 
20,  21,  22,  24,  25,  28,  31,  35,  36,  37.  Therefore,  (3.8)  leads  to  99  affine  inequality  constraints.  This 
problem  is  readily  solvable  by  1  inprog  in  MATLAB’s  Optimization  Toolbox  [11].  The  active 
constraints  are  checked  with  a  MATLAB  script  to  determine  the  optimal  control  policy. 

Figure  3.1  shows  an  example  of  4  stationary  control  policies  depicted  in  terms  of  indicator 
functions,  as  defined  in  (3.4),  of  the  state  x  G  X.  It  can  be  seen  that  the  optimal  policy  takes  into 
consideration  of  anticipated  targets  more  than  the  aggressive  policy,  but  is  much  more  aggressive 
in  terms  of  use  of  resources  than  the  conservative  policy.  Among  the  four  optimal  policies  solved, 
only  the  policy  derived  under  the  least  sensor  loss  is  plotted  (bottom),  which  will  be  shown  shortly 
to  outperform  other  three  optimal  policies  in  terms  of  both  MTTNE  and  availability.  The  minimum 
queue  length  policy  coincides  with  the  greedy  policy,  as  expected.  The  other  two  optimal  policies 
make  less  aggressive  use  of  resources  than  the  optimal  policy  shown.  The  least  sensor  loss  policy 
will  be  called  the  optimal  policy  from  this  point  on. 

The  control  policies  are  robust  with  respect  to  the  range  of  parameter  variations  that  have  been 
examined:  zg  G  [0.001,0.01]  1/sec.,  and  A  G  [0.001,0.1]  1/sec.  The  optimal  policy  is  calculated 
at  a  =  0.0792.  All  optimal  policies  drift  slightly  toward  more  conservative  actions  (using  fewer 
resources)  as  the  discount  factor  a  increases,  which  is  consistent  with  the  outcome  of  the  longer 
term  policy  making.  Because  of  the  finite  target  population  setup,  the  effect  of  increasing  the  target 
arrival  rate  is  not  fully  reflective  of  the  target  traffic  intensity.  Simulations  with  MATLAB  SimEvets 
[10]  are  being  performed  without  limiting  the  target  population. 
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Figure  3.1:  Control  policy  indicators,  red:  one  sever/target;  green:  two  servers/target;  blue:  three 
servers/target 
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4  Performance  analysis 


The  seven  policies  developed  in  Section  3  are  compared  against  one  another  with  respect  to  two 
common  measures  of  fault-tolerance:  mean  time  to  network  expiration  (MMTNE)  and  availability. 
These  have  been  used  in  [15]  in  a  similar  fashion  as  performance  measures  of  a  database  unit. 

When  no  replenishment  is  provided,  the  network  life  eventually  expires  when  all  sensor-pairs 
are  lost.  This  occurs  when  the  network  enters  one  of  its  absorbing  states  at  34,  35,  36,  or  37. 
Decompose  the  state  probability  vector 

7 r(f)  =  7TqC0]  (4.1) 

1x33  1x4 

where  vector  7Tr(f)  contains  transient  state  probabilities,  and  n a(t)  contains  absorbing  state  proba¬ 
bilities.  Decomposing  the  rate  transition  matrix  Q  accordingly  yields 


(  )  _  Q ii  Q 12 
[  0  0 


(4.2) 


From  (4.2),  it  can  be  determined  that  mean  time  to  network  expiration  is  given  by 

MTTNE  =  -vrT(0)g711lT,  1T  =  [1__^_JJT  (4.3) 

1x33 

Suppose  as  soon  as  the  network  expires,  a  replenishment  process  starts.  Suppose  with  a  rate 
ijj  the  airborne  sensors  are  replenished,  and  at  the  completion  of  the  replenishment,  the  tasking 
process  immediately  resumes.  In  this  case,  the  Markov  chain  (2.5)  becomes  irreducible,  and  a 
unique  steady-state  distribution  exists  [9].  The  steady-state  availability,  which  can  be  roughly 
thought  of  as  the  fraction  of  time  the  network  has  at  least  one  surviving  pair  of  sensors,  is  computed 
by 

Anet  =  1  -7TF(oo),  (4.4) 

where  ttf(oo)  =  7r34(oo)  +  7r35(oo)  +  7r36(oo)  +  7r37(oo),  the  sum  of  state  probabilities  associated 
with  network  expiration,  which  can  be  determined  by  solving 

37 

n(oo)Q  =  0,  and  ^7rx(oo)  =  1.  (4.5) 

X=1 
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A  slightly  different  notion  of  availability  is  also  examined,  where  the  network  is  considered 
unavailable  as  long  as  unattended  targets  are  present.  In  this  case,  the  network  availability  is  given 

by 


Atgt  ~ 

11 

^7Ti(00)  +  7Ti3(oo)  +  7Ti6(oo)  +  7T2(oo) 
i=l 

+7r23(oo)  +  7T24(oo)  +  7T26(oo)  +  vr27(oo).  (4.6) 


co=0,  u  =1/60,  p  =1 .5p  ,  p  =1 .5p  ,  v  =0.0005,  v  =0.001  to  0.01  ,v  =1 ,5v  ,v  =1 .5vr 
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Figure  4.1:  Mean  time  to  network  expiration  as  a  function  of  target  arrival  rate  with  two  sets  of 
sensor  loss  rates  as  parameters. 

Mean  time  to  network  expiration  is  plotted  in  Figure  4.1  against  target  arrival  rate  with  two  sets 
of  sensor  loss  rates  as  parameters  at  a  =  0.0792.  It  shows  that  compromise  that  optimal  policy 


12 


makes  between  being  too  greedy  and  too  conservative  enhances  the  network  life  consistently  for 
all  parameters  values  considered. 
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0.01,  v2=1.5v1,v3=1.5v2 


0.5 
0.48 
0.46 
0.44 
0.42 
I  0.4 
0.38 
0.36 
0.34 
0.32 


3 


Figure  4.2:  Network  availability  with  at  least  one  surviving  server. 


In  Figure  4.2,  availability  is  also  plotted  against  target  arrival  rate  with  two  sets  of  sensor  loss 
rates  as  parameters  at  discount  factor  a  =  0.0792.  The  observations  from  the  MTTNE  apply  in 
terms  of  the  gain  the  optimal  policy  offers.  It  is  noted  that  the  availability  is  low.  This  is  because 
of  the  low  replenishment  rate  used  in  the  computation,  which  corresponds  to  an  expected  time  of 
more  than  40  minutes  to  reestablish  the  lost  network. 

It  is  expected  that  the  network  availability  defined  as  the  probability  that  all  targets  is  lower 
than  the  availability  defined  as  the  probability  that  there  is  at  least  one  surviving  server.  On  the 
other  hand,  the  dependence  of  both  notions  of  availability  on  the  sensor  loss  rate  and  on  the  target 
arrival  rate  stays  the  same.  Figure  4.3  shows  the  plot  of  the  two  availabilities  against  target  arrival 
rate  with  two  sets  of  sensor  loss  rates  ul  as  parameters  under  the  optimal  policy. 
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Extensive  simulations  using  [10]  is  being  conducted  to  generate  a  more  complete  picture  of  the 
network  performance  in  response  to  control  policies  for  a  larger  size  of  networks.  The  results  will 
be  reported  elsewhere. 
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Availiability 


co=0.0004,  jll1 =1/60,  p.  =1.5^,  jli3=1  -5jli2,  v0=0.0005, 
v^O.OOI  to  0.01 ,  v2=1 .5v1 ,  v3=1  ,5v2 
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Figure  4.3:  Comparison  between  availability  with  at  least  one  surviving  server  (dash)  and  avail¬ 
ability  with  all  targets  being  attended  (solid). 
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5  Conclusions 


This  work  sought  to  determine  supervisory  control  policies  that  best  configure  an  airborne  location 
sensor  network  to  provide  a  high  degree  of  guarantee  of  prompt  completion  of  coordinated  data 
acquisition  and  processing  missions  in  the  face  of  loss  of  vehicles.  Use  of  redundancy  and  dynamic 
allocation  of  participating  sensors  was  the  key  enablers. 

The  report  presented  a  queuing  network  approach  to  optimal  sensor-pair  assignment  to  locate 
detected  targets  for  a  small  scale  airborne  sensor  network,  where  use  of  redundancy  is  balanced 
with  avoiding  more  vehicle  exposure. 

A  number  of  related  issues  are  being  investigated.  In  addition  to  loss  of  sensors,  degradation 
of  network  performance  can  also  be  the  result  of  broken  communication  links.  Such  incidents  are 
modeled  as  intermittent  faults  of  the  servers  in  queuing  networks.  The  effects  of  such  faults  will  be 
studied  using  discrete  event  simulations,  which  will  also  examine  possible  emergent  phenomenon 
of  larger  networks,  and  non-homogeneous  sensors  and  targets. 

A  highly  relevant  task  is  to  solve  a  guidance  and  then  a  control  problem  of  the  vehicles.  Guid¬ 
ance  problem  [7]  refers  to  that  of  establishing  a  criterion  and  deriving  a  set  of  desired  vehicle 
trajectories  under  the  criterion  that  the  airborne  sensors  are  expected  to  follow  to  expedite  the  tar¬ 
get  location  estimation  to  a  required  accuracy.  It  is  known  that  the  quality  of  acquired  data  by 
the  airborne  sensors  depends  highly  on  both  the  network  architecture  which  is  determined  by  the 
number  of  sensor-pairs,  and  the  states  of  all  participating  sensors  relative  to  the  target  and  to  one 
another.  A  guidance  principle  based  on  the  entropy  [5]  of  the  noise  distribution  of  the  sensed  sig¬ 
nal  has  been  established,  based  on  which  one  seeks  to  adjust  the  states  of  the  sensors  to  the  most 
suitable  positions  and  velocities  for  further  data  acquisition  and  processing. 

Once  the  guidance  principle  is  determined,  feasible  vehicle  trajectories  can  be  generated.  Path 
following  control  can  be  performed  with  time  coordination  [8]  to  achieve  synchronous  data  acqui¬ 
sition  by  the  sensors.  Both  the  guidance  and  the  control  problems  are  being  investigated,  and  will 
be  reported  separately  in  the  near  future. 
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6  Designs  and  simulation  modeling 


In  general,  the  sensors  of  the  sensor  networks  are  subject  to  random  breakdowns.  This  has  a  heavy 
influence  on  the  performance  measures.  Thus,  if  we  model  a  system  containing  unreliable  sensor- 
pairs,  it  is  important  to  take  it  into  account  in  the  model  construction.  Of  course,  the  breakdown 
of  the  sensor-pairs  has  the  most  significant  negative  impact  on  the  performance  of  the  system.  The 
complexity  of  the  system,  such  as  its  stochastic,  unpredictable  behavior,  independencies  between 
individual  events  and  the  informal  nature  of  many  events,  makes  the  analysis  profoundly  problem¬ 
atic.  Such  a  modeling  problem  is  very  significant  in  large  complex  systems.  The  advantage  of 
using  simulation  is  listed  below: 

1.  Simulation  provides  a  controlled  environment  for  performance  measurements. 

2.  Changes  can  be  made  easily  to  the  model. 

3.  Simulation  and  modeling  provides  an  easy  way  to  control  influence  of  parameters  and  con¬ 
duct  sensitivity  analysis. 

4.  It  does  not  rely  on  mathematical  concepts  heavily  and  is  a  very  robust  technique. 

Simulation  of  the  airborne  system  is  conducted  in  Simulink  of  MATLAB  using  the  event-based 
package  Simevents  [10],  [11].  We  began  with  three  simulation  models.  These  three  models  are 
closed  queue  with  finite  homogenous  sensor-pairs  and  finite  arrival  targets.  The  first  model  is 
reducible,  without  both  system  replenishment  and  channel  fading.  We  used  this  model  to  measure 
the  mean  time  to  network  expiration  (MTTNE).  The  second  model  is  irreducible  with  the  feature 
of  system  replenishment.  The  third  model  is  also  irreducible  with  both  replenishment  and  channel 
fading.  Both  the  second  and  third  model  is  used  to  measure  steady  state  availability  and  mean 
response  time  of  the  system.  Several  control  policies  associated  with  tasking  are  examined  in  the 
simulation  and  compared  with  the  analytic  results. 


6.1  Basic  algorithm  of  the  simulation  models 

Events  do  not  have  equal  probability;  some  events  are  more  likely  to  occur  than  others;  the  exe¬ 
cution  order  of  simultaneous  events  is  set  to  be  randomized  instead  of  arbitrary  through  all  simu¬ 
lations.  Upon  the  occurrence  of  permanent  failures,  targets  are  ejected  out  from  the  sensor-pairs 
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before  service  completion  and  re-attending  the  queue  with  priority,  and  they  require  the  same  ser¬ 
vice  time  that  was  assigned.  These  targets  preempt  the  waiting  targets  in  the  queue  but  don’t 
preempt  the  targets  that  are  already  in  service.  Upon  occurrence  of  intermittent  failures,  targets  are 
ejected  out  from  the  sensor-pairs,  and  they  memorize  their  residual  event  lives.  They  immediately 
return  to  the  assigned  sensor-pairs  upon  recovery  of  the  sensor-pairs.  The  system  is  said  to  be 
down  when  all  sensor-pairs  end  their  life.  No  actions  can  be  taken  during  the  repairing  process. 
The  event  of  failure  is  activated  again  once  the  system  is  operational.  The  system  contains  a  finite 
population  of  targets;  the  delay  elements  in  the  feedback  loop  generate  new  targets  upon  service 
completion  of  the  targets  from  one  or  more  sensor-pairs  from  the  upper  part  of  the  closed  loop 
system.  Appendix  C  convolutes  the  simulation  and  implements  the  following  algorithm: 

1.  Initialize  the  system  to  state  0000. 

2.  Activate  all  feasible  events 

3.  Create  an  initial  entity  and  set  a  service  time  for  it. 

4.  Service  time  is  updated  according  to  the  number  of  entities  in  the  system. 

5.  Failure  time  of  the  servers  is  updated  with  respect  to  the  number  of  entities  in  the  system  as 
well.  A  new  failure  entity  will  be  generated  by  the  random  number  generator  and  preempts 
the  previous  failure  entity  when  the  queue  length  changes.  The  old  failure  entity  is  discarded 
under  this  condition. 

6.  Compare  the  transition  times  of  the  feasible  events;  choose  the  shortest  time  and  resolve  the 
trigger  event.  Update  to  new  state  upon  service  completion,  entity  arrival  and  server  failure. 

7.  Simulation  stops  for  non-repairable  systems  once  the  absorbing  states  are  reached.  For  re¬ 
pairable  systems,  simulation  stops  only  if  there  is  a  stopping  rule. 

For  most  communication  systems,  redundancy  is  desired.  Assume  the  stopping  rule  of  this 
airborne  system  is  defined  as  the  system  is  down  when  the  backups  of  severs  are  also  gone  or  when 
the  system  is  only  allowed  to  replenish  once.  Another  issue  has  to  do  with  the  failure  time.  If  the 
current  simulation  time  is  greater  than  the  failure  time  of  a  sensor-pair,  the  sensor-pair  fails  right 
away.  If  the  current  simulation  time  is  less  than  the  failure  time,  then  the  service  time  of  the  target 
is  equal  to  the  failure  time  minus  the  old  simulation  time. 


6.2  Techniques  for  steady  state  analysis 

Unlike  in  queuing  theory  where  steady  state  results  in  analytic  models  are  easily  obtainable,  the 
steady  state  simulation  is  not  an  easy  task.  Gathering  steady  state  simulation  output  requires  a 
statistical  assertion  that  the  simulation  model  reached  the  steady  state.  The  main  difficulty  is  to 
obtain  independent  simulation  runs  with  exclusion  of  the  transient  period.  There  are  two  techniques 
commonly  used  for  steady  state  simulation,  the  Method  of  Batch  means,  and  the  Independent 
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Replication  [18],  [19].  None  of  these  two  methods  is  superior  to  the  other  in  all  cases.  Intuitively 
one  may  say  the  method  of  independent  replication  is  superior  in  producing  statistically  good 
estimates  for  the  system’s  performance  measurements.  In  fact,  not  one  method  is  superior  in  all 
cases,  and  it  all  depends  on  the  traffic  intensity. 

We  are  using  independent  replication  in  this  experiment.  This  technique  gets  n  independent 
runs  of  the  simulation  experiment  by  running  the  simulation  n  times  with  different  initial  random 
seeds  for  the  simulator’s  random  number  generator.  The  main  question  is  determination  of  the 
number  of  runs.  The  confidence  level  of  simulation  output  drawn  from  a  set  of  simulation  runs 
depends  on  the  size  of  the  data  set.  The  larger  the  number  of  runs,  the  higher  is  the  associated 
confidence.  However,  more  simulation  runs  require  more  effort  for  large  systems.  Hence,  the 
main  goal  is  to  find  the  smallest  number  of  simulation  runs  that  will  also  provide  the  desired  con¬ 
fidence.  One  main  limitation  of  these  designs  is  that  the  outputs  are  random  from  the  simulations. 
Thus,  estimates  of  parameter  effects  are  subject  to  possibly  considerable  variance.  Unlike  phys¬ 
ical  experiments,  we  have  the  luxury  in  simulation  of  replicating  the  runs  many  times  to  reduce 
this  variance  or  perhaps  replicating  the  whole  design  many  times  to  get  many  independent  and 
identically  distributed  estimates  of  main  parameter  effects,  which  could  then  combine  to  form  a 
confidence  interval.  Due  to  the  lengthy  simulation  time,  all  results  are  averaged  over  30  simulation 
runs. 


6.3  Performance  analysis  from  simulation 

To  run  a  simulation  multiple  times  and  gather  statistics,  the  initial  seed  value  would  have  to  be 
reassigned  for  each  run  in  all  blocks  containing  it.  The  MATLAB  code  initial  seed  value  generator 
from  [10]  shown  in  appendix  9.4  generates  random  initial  seeds  for  the  system.  This  is  one  of  the 
methods  that  may  be  used  to  obtain  results  from  multiple  runs.  To  see  the  performance  over  the 
range  of  arrival  rate,  we  can  also  use  the  code  ’Running  a  simulation  and  varying  a  parameter’ 
from  appendix  9.4.  The  following  subsections  are  the  performance  analysis  with  data  collected 
from  simulations.  They  are  plotted  in  Excel. 

6.3.1  Closed  -  queue  reducible  model 

Fig  6.1  shows  the  Mean  time  to  network  expiration  vs.  arrival  rate.  It  has  the  same  arrangement  as 
fig.  3.3  with  the  closeness  of  the  curves  at  slow  arrival  rate  and  the  spread  of  the  curves  at  faster 
arrival  rate.  The  result  is  not  exactly  the  same  as  the  analytic  result,  but  very  close. 
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Figure  6.1:  Mean  time  to  network  expiration  vs.  arrival  rate 


(simulation  model  1) 
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6.3.2  Closed  -  queue  irreducible  model 


Figure  6.2:  Expected  response  time  vs.  arrival  rate  (simulation  model  2) 
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Figure  6.3:  Availability  vs.  arrival  rate  (simulation  model  2) 

Compare  fig.  6.3  to  fig  3.4,  the  system  availability  obtain  from  the  simulation  is  much  higher 
than  the  analytic  model,  that  is  because  the  simulation  only  allows  to  replenish  the  system  once. 
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6.3.3  Closed  -  queue  irreducible  model  with  intermittent  channel  fading 


Figure  6.4:  Expected  response  time  vs.  arrival  rate  (simulation  model  3) 

From  fig.  6.4,  the  response  time  for  the  system  with  intermittent  channel  fading  is  longer  than 
the  response  time  without  intermittent  channel  fading.  It  is  expected  that  the  time  to  locate  a  target 
would  be  extended. 


6.4  Simulation  conclusion 

The  set  up  of  the  simulation  models  is  somewhat  different  from  the  analytic  model.  First  of  all, 
in  the  analytic  model,  the  performance  measurements  are  obtained  under  all  possible  states  of  the 
system.  In  the  simulation  models,  the  simulation  is  running  under  the  selection  of  the  control 
policies,  not  all  system  states  would  be  presented.  For  example,  if  the  conservative  policy  is 
applied,  a  state  in  which  all  three  sensor-pairs  tie  together  to  serve  a  target  would  never  occur. 

Secondly,  change  of  simulation  results  may  be  due  to  the  run  time  of  a  model  and  randomness 
of  the  event  lives.  Compare  fig.  6.1  to  3.3,  the  MTTNE  is  slightly  different  from  the  Markov 
model.  The  expected  response  times  obtained  from  simulation  model  3  are  longer  than  those  for 
model  2.  This  is  predictable  where  the  intermittent  channel  fading  extends  the  time  to  locate  targets 
(additional  time  is  needed  for  data  recovery).  The  response  times  from  the  simulation  models  2 
and  3  are  not  the  same  as  fig.  3.6b  because  fig  3.6b  is  plotted  under  the  assumption  that  the  number 
of  served  targets  for  all  policies  is  constant.  The  number  of  served  targets  in  simulation  models 
is  correct  and  changes  by  the  arrival  and  service  rate.  Also,  another  fact  is  that  the  run  time  for 
simulations  was  too  short;  the  systems  are  not  under  steady  state  condition.  A  more  ideal  set  of 
data  can  be  obtained  by  having  longer  run  time  with  a  larger  number  of  departed  targets  over  time. 
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Figure  6.5:  Availability  vs.  arrival  rate  (simulation  model  3) 


Similarly,  the  availabilities  from  the  simulation  models  are  higher  than  the  analytic  model.  This  is 
also  due  to  the  short  simulation  run  time  and  the  effect  of  non-steady  state  operation. 
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8  MATLAB  code  for  the  analytic  model  of 
the  airborne  network 


8.1  Transition  rate  matrix 

The  transition  rate  matrix  is  computed  from  the  transition  rate  diagram  from  the  previous  page  or 
from  Table  4. 


8.2  Mean  time  to  network  expiration 


clc,  clear 
%  parameters 

la=linspace  (1/300,  1/50,  30); 
nul  =  linspace  (1/1000,  1/100,  2); 
nu0=0.0005;  nu2=1.5*nul;  nu3=1.5*nu2; 
mul=l/60;  mu2=1.5*mul;  mu3=1.5*mu2; 

%%  compute  mean  time  to  network  expiration  (reducible  case) 
for  i=l : length (la) 

for  j=l : length (nul ) 

setpara(mul,  mu2,  mu3,  nuO,  nul(j),  nu2(j),  nu3(j),  la(i),  0); 

Q  =Qnewl; 

pi_T  =  [1  zeros (1,  (length (Q) -5 ) )  ] ; 
ones_T  =  ones  (  (length (Q) -4 ), 1 )  ; 

Q_T  =  Q (1  :  (length (Q) -4) , 1 :  (length (Q) -4) ) ; 

MTTNEl(i,j)  =  -pi_T  *  inv(Q_T)  *  ones_T;  %  greedy 

Q2=Qnew2 ; 

pi_T2  =  [1  zeros ( 1 ,( length (Q2 ) -5 ))] ; 
ones_T2  =  ones ( (length (Q2 ) -4 ), 1 ) ; 

Q2_T  =  Q2 (1 : (length (Q2) -4) , 1 : (length (Q2) -4) ) ; 

MTTNE2(i,j)  =  -pi_T2  *  inv(Q2_T)  *  ones_T2;  %  conservative 
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Qag=Qnewag; 


pi_Tag  =  [1  zeros (1,  ( length (Qag) -5 ))] ; 
ones_Tag  =  ones ( (length (Qag) -4 ) , 1 ) ; 

Qag_T  =  Qag (1 : (length (Qag) -4) , 1 : (length (Qag) -4) ) ; 

MTTNEag(i, j)  =  -pi_Tag  *  inv(Qag_T)  *  ones_Tag;  %  aggressive 


Qopt l=Qnewopt 1 ; 

pi_Toptl  =  [1  zeros (1,  (length (Qoptl) -5) )] ; 
ones_Toptl  =  ones ( (length (Qoptl) -4) , 1 )  ; 

Qoptl_T  =  Qoptl (1 : (length (Qoptl) -4) , 1 : (length (Qoptl) -4) ) ; 
MTTNEopt 1 (i, j )  =  -pi_Toptl  *  inv(Qoptl_T) 

*  ones_Toptl;  %  optimal  policy  1 


Qopt2=Qnewopt2 ; 

pi_Topt2  =  [1  zeros (1,  (length (Qopt2) -5) )] ; 
ones_Topt2  =  ones ( (length (Qopt2 ) -4 ) , 1 ) ; 

Qopt2_T  =  Qopt2 (1 :  (length (Qopt2) -4) , 1 :  (length (Qopt2) -4) ) ; 
MTTNEopt 2 (i, j)  =  -pi_Topt2  *  inv(Qopt2_T) . . . 

*  ones_Topt2;  %  optimal  policy  2 


Qopt  3=Qnewopt 3 ; 

pi_Topt3  =  [1  zeros (1,  (length (Qopt3) -5) )] ; 
ones_Topt3  =  ones  (  (length (Qopt3 ) -4 ) , 1 ) ; 

Qopt3_T  =  Qoptl (1 : (length (Qopt3) -4) , 1 : (length (Qopt3) -4) ) ; 
MTTNEopt 3 (i,  j )  =  -pi_Topt3  *  inv (Qopt 3_T)  .  .  . 

*  ones_Topt3;  %  optimal  policy  3 

end, end 


diffl  =  MTTNEag-MTTNEl; 
dif f 2  =  MTTNE1-MTTNE2; 
dif f 3  =  MTTNEag  -  MTTNE2 ; 


%  diff  between  aggressive  &  greedy 
%  diff  between  greedy  &  conservative 
%  diff  between  aggressive  &  conservative 


%%  plot  MTTNE  vs .  arrival  rate 


figure ( 1 ) 
subplot (1,2,1) 

plot ( la, MTTNE2 ( : , 1) ,  la, MTTNE1 ( : , 1) ,  ’ -x ' ,  la, MTTNEag ( : , 1) ,  la, MTTNEopt 1 ( : , 1) ) 

title ({ 1 \omega=0,  \mu_l=l/60,  \mu_2=l . 5\mu_l ,  \mu_3=l . 5 \mu_2 ,  \nu_0=0 . 0 0 05 , . . . 
\nu_l  =  0.001  to  0.01,  \nu_2=l . 5\nu_l ,  \nu_3=l . 5\nu_2 '  ;  .  .  . 

’MEAN  TIME  TO  NETWORK  EXPIRATION  VS.  ARRIVAL  RATE ' ; ' \nu_l  -  0.001’}) 


legend ( ' conservative ' , ' greedy ' , ’ aggressive ' , ' optimal  policy  1 ' ) 
xlabel (' Arrival  Rate  (\lambda  1/sec.) ') 
ylabel ( 'MTTNE  (sec.)'),  xlim( [1/300,  1/50]) 
subplot (1,2,2) 

plot ( la, MTTNE2 ( : , 2) ,  la, MTTNE1 ( : , 2) ,  ' -x ' ,  la, MTTNEag ( : , 2) ,  la, MTTNEopt 1 ( : , 2) ) 

title ( ' \nu_l  =  0.01'),  xlabel (' Arrival  Rate  (\lambda  1/sec. ) ' ) 
ylabel ( 'MTTNE  (sec.)'),  xlim( [1/300,  1/50]) 
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figure  (2 ) 

subplot  (1,2, 1) ,  plot (la, MTTNEopt 1  (  : , 1 ) ,  '-+',  la,  MTTNEopt2  (  : , 1 ) ,  '-o’,.... 

la, MTTNEopt3  (  : , 1) ) 

title ({' \omega=0,  \mu_l=l/60,  \mu_2=l . 5\mu_l ,  \mu_3=l . 5\mu_2 , . . . 

\nu_0=0 . 0005,  \nu_l=0.001  to  0.01,  \nu_2=l . 5\nu_l ,  \nu_3=l . 5\nu_2 ' ; . . . 
’MEAN  TIME  TO  NETWORK  EXPIRATION  VS.  ARRIVAL  RATE ’ ; ’ \nu_l  =  0.001’}) 

legend (' optimal  policy  1',  'optimal  policy  2',' optimal  policy  3') 
xlabel (' Arrival  Rate  (\lambda  1/sec.)'),  ylabel ( ' MTTNE  (sec.)') 
xlim ( [1/300,  1/50]  ) 
subplot (1,2,2) 

plot (la, MTTNEopt 1 ( : , 2 ) ,  la,  MTTNEopt2 ( : , 2 ) ,  '-o',  la, MTTNEopt 3 ( : , 2 ) ) 

title  (' \nu_l  =  0.01') 

xlabel (' Arrival  Rate  (\lambda  1/sec.) ') 
ylabel ( 'MTTNE  (sec. )  '  ) 
xlim ( [1/300,  1/50] ) 

8.3  Steady-state  availability 


clc, clear 


%  parameters 

la=linspace  (1/300,  .02,  30); 

nul  =  linspace  (1/1000,  1/100,  2); 
nu0=0.0005;  nu2=1.5*nul;  nu3=1.5*nu2; 
mul=l/60;  mu2=1.5*mul;  mu3=1.5*mu2; 
w=0 . 0004 ; 


%%  computing  steady  state  availability 
for  i=l : length ( la) 

for  j=l : length (nul ) 

setpara(mul,  mu2,  mu3,  nuO,  nul(j),  nu2(j),  nu3(j),  la(i),  w) ; 


Q  =Qnewl; 
Q2=Qnew2 ; 
Qag=Qnewag; 
Qoptl=Qnewoptl ; 
Qopt2=Qnewopt2 ; 
Qopt3=Qnewopt3; 


%greedy 
% conservative 
%aggressive 
%optimal  1 
%optimal  2 
%optimal  3 


Q (:, length (Q) )  =  [ones ( 1 , length (Q) )] ; 
p_ss  =  [ zeros ( 1 ,( length (Q) -1 )) ,  l]*inv(Q); 


%avail  1  when  the  system  is  up 

A_sysl(i,j)  =  p_ss (1 : (length (Q) -4) ) *ones ( (length (Q) -4) , 1) ; 


%avail  2:  all  targets  are  in  service 

A_sysl_2 ( i, j )  =  l-p_ss ( 12 ) -p_ss ( 14 ) -p_ss ( 15 ) -p_ss ( 17 ) -p_ss ( 18 ) -p_ss (19) . . . 
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-p_s  s ( 2  0 ) -p_s s ( 2 1 ) — p s s ( 2  5 ) -p_s  s ( 2  8 ) -p_s  s ( 2  9 ) -p_s  s ( 3  0 ) -p_s  s ( 3 1 )  . . . 

-p_s  s  ( 3  2 ) -p_s  s ( 3  3 ) — p _ s s ( 3  4 ) -p_s s ( 3  5 ) -p_s s ( 3  6 ) -p_s s (37) ; 

Q2 (:, length (Q2) )  =  [ones ( 1 , length (Q2 ))] ; 
p_ss2  =  [zeros (1,  (length (Q2) -1)  )  ,  l]*inv(Q2); 

A_sys2(i,j)  =  p_ss2 ( 1 : (length (Q2 ) -4 )) *ones ( (length (Q2 ) -4 ), 1 ) ; 

A_sys2_2 (i, j)  =  l-p_ss2  (12) -p_ss2 (14) -p_ss2 (15) -p_ss2 (17) -p_ss2 (18)  . . . 

-p_s s  2 (19) -p_s s  2 (20) -p_s  s  2 (21) -p_s  s  2 (25) -p_s  s  2 (28) -p_s  s  2 ( 2  9 )  . .  . 

-p_s  s  2 (30) -p_s s  2 (31) -p_s  s  2 (32) -p_s  s  2 (33) -p_s  s  2 (34) -p_s  s  2 ( 3  5 )  . .  . 

-p_s s  2 (36) -p_s  s  2 (37) ; 

Qag (: , length (Qag) )  =  [ones  (1, length (Qag) )] ; 
p_ssag  =  [zeros (1,  (length (Qag) -1) ) ,  l]*inv(Qag); 

A_sysag (i,  j)  =  p_ssag (1 :  (length (Qag) -4) ) *ones ( (length (Qag) -4) , 1) ; 
A_sysag_2 ( i , j )  =  l-p_ssag (12) -p_ssag (14) -p_ssag (15) -p_ssag (17) -p_ssag (18) . . . 

-p_ssag (19) -p_ssag (20) -p_ssag (21) -p_ssag (25) -p_ssag (28) -p_ssag (29) . . . 
-p_ssag (30) -p_ssag (31) -p_ssag (32) -p_ssag (33) -p_ssag (34) -p_ssag (35) . . . 
-p_s  s  ag ( 3  6 ) -p_s  s  ag ( 3  7 )  ; 

Qoptl (:, length (Qoptl ) )  =  [ ones ( 1 , length (Qoptl ))] ; 
p_ssoptl  =  [zeros (1, (length (Qoptl) -1) ) ,  1 ] *inv (Qoptl) ; 

A_sysopt 1 (i, j )  =  p_ssopt 1 ( 1 : (length (Qoptl ) -4 )).. . 

★ones ( (length (Qoptl ) -4) , 1) ; 

A_sysoptl_2 (i, j )  =  l-p_ssoptl (12) -p_ssoptl (14) -p_ssoptl (15) -p_ssoptl (17) . . . 

-p_ssoptl (18) -p_ssoptl (19) -p_ssoptl (20) -p_ssoptl (21) -p_ssoptl (25) . . . 
-p_ssopt 1 (28 ) -p_ssoptl (29) -p_ssopt 1 (30 ) -p_ssoptl (31 ) -p_ssoptl (32 ) . . . 
-p_ssoptl (33) -p_ssoptl (34) -p_ssoptl (35) -p_ssoptl (36) -p_ssoptl (37) ; 

Qopt2 (:, length (Qopt2 ) )  =  [ ones ( 1 , length (Qopt2 ))] ; 
p_ssopt2  =  [zeros (1,  (length (Qopt2) -1) ) ,  1 ] *inv (Qopt2 ) ; 

A_sysopt2 (i, j )  =  p_ssopt2 ( 1 :( length (Qopt2 ) -4 )).. . 

★  ones ( (length (Qopt2 ) -4 )  ,  1)  ; 

A_sysopt2_2 (i, j )  =  l-p_ssopt2 (12) -p_ssopt2 (14) -p_ssopt2 (15) -p_ssopt2 (17) . . . 

-p_ssopt2 (18) -p_ssopt2 (19) -p_ssopt2 (20) -p_ssopt2 (21) -p_ssopt2 (25) . . . 
-p_ssopt2  (28 ) -p_ssopt2 (2  9) -p_ssopt2 (30 ) -p_ssopt2 (31 ) -p_ssopt2 (32 )  .  .  . 
-p_ssopt2 (33) -p_ssopt2 (34) -p_ssopt2 (35) -p_ssopt2 (36) -p_ssopt2 (37) ; 


Qopt 3 (:, length (Qopt3 ) )  =  [ ones ( 1 , length (Qopt3 ))] ; 
p_ssopt3  =  [zeros  (1,  (length (Qopt3) -1) ) ,  1 ] *inv (Qopt3) ; 

A_sysopt3 (i, j )  =  p_ssopt3 ( 1 : ( length (Qopt 3 ) -4 ) ) . . . 

★ones ( (length (Qopt 3) -4) , 1) ; 

A_sysopt3_2 (i, j )  =  l-p_ssopt3 (12) -p_ssopt3 (14) -p_ssopt3 (15) -p_ssopt3 (17) . . . 

-p_ssopt3 (18) -p_ssopt3 (19) -p_ssopt3 (20) -p_ssopt3 (21) -p_ssopt3 (25) . . . 
-p_ssopt3 (28) -p_ssopt3 (29) -p_ssopt3 (30) -p_ssopt3 (31) -p_ssopt3 (32) . . . 
-p_ssopt3 (33) -p_ssopt3 (34) -p_ssopt3 (35) -p_ssopt3 (36) -p_ssopt3 (37) ; 

end, end 

dif f l=A_sysag  -  A_sysl;  %  diff  between  aggressive  and  greedy 

diff2=A_sysl  -  A_sys2;  %  diff  between  greedy  and  conservative 
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dif f 3=A_sysag-A_sys2 ;  %  diff  between  aggressive  and  conservative 

%%  plots:  avail  1  vs.  arrival  rate 

figure ( 1 ) 
subplot (1,2,1) 

plot (la, A_sys2  (  : , 1) ,  la, A_sysl  (  : , 1 ) ,  ' -x ' ,  la, A_sysag  ( : , 1 ) ,  la,  A_sysopt 1  (  : , 1 ) ) 
title ({' \omega=0 . 0004 ,  \mu_l  =  l/60,  \mu_2  =  l . 5 \mu_l ,  \mu_3=l . 5\mu_2 ,  . .  . 

\nu_0  =  0 . 0005,  \nu_l=0.001  to  0.01,  \nu_2  =  l . 5\nu_l ,  \nu_3=l . 5\nu_2  ' ;  . . . 
'AVAILABILITY  VS.  ARRIVAL  RATE  1 ;  ' \nu_l  -  0.001'}) 
legend ( ' conservative ' , ' greedy ' , ' aggressive ' , ' optimal  policy  1 ' ) 
xlabel (' Arrival  Rate  (\lambda  1 /sec . ) ' ) , ylabel ( ' Availiability ' ) 
xlim ( [1/300,  1/50] ) 
subplot (1,2,2) 

plot (la, A_sys2  (  : , 2) ,  la, A_sysl  (  : , 2 ) ,  ' -x ' ,  la, A_sysag  ( : , 2 ) ,  la,  A_sysopt 1  (  : , 2 ) ) 
title (' \nu_l  =  0.01'),  xlabel (' Arrival  Rate  (\lambda  1/sec.)') 
ylabel (' Availiability ') ,  xlim( [1/300,  1/50]) 

figure  (2 ) 

subplot ( 1 , 2 , 1 ) ,  plot (la, A_sysoptl ( : , 1 ) ,  la, A_sysopt2 ( : , 1 ) ,  ' -+ ' ,  la, A_sysopt3 ( : , 1 ) ) 

title ({' \omega=0 . 0004 ,  \mu_l=l/60,  \mu_2=l . 5 \mu_l ,  \mu_3=l . 5\mu_2 , . . . 

\nu_0  =  0 . 0005 ,  \nu_l  =  0.001  to  0.01,  \nu_2  =  l . 5\nu_l ,  \nu_3=l . 5\nu_2 ' ;  . . . 
'AVAILABILITY  VS.  ARRIVAL  RATE ' ; ' \nu_l  =  0.001'}) 
legend (' opt imal  policy  1',  'optimal  policy  2', 'optimal  policy  3') 
xlabel (' Arrival  Rate  (\lambda  1/sec.)'),  ylabel (' Availiability ' ) 
xlim ( [1/300,  1/50]  ) 

subplot ( 1 , 2 , 2 ) ,  plot (la, A_sysoptl ( : , 2) ,  la, A_sysopt2 ( : , 2 ) ,  ' -+ ' ,  la, A_sysopt3 ( : , 2 ) ) 

title ( ' \nu_l  =  0.01'),  xlabel (' Arrival  Rate  (\lambda  1/sec.)') 
ylabel (' Availiability ')  ,  xlim( [1/300,  1/50]) 

%%  plot  of  avail  1  vs.  avail  2 
figure  ( 3 ) 
subplot  (1,2, 1 ) 

plot (la, A_sysopt 1 ( : , 1 ) , ' - . r ' ,  la, A_sysopt2 ( : , 1 ) ,  ' -xg ' , . . . 

la, A_sysopt3 ( : , 1 ) , ' -ob ' ,  la, A_sysopt 1_2 ( : , 1 ) , ' r ' ,  la,... 

A_sysopt2_2 ( : , 1 ) ,  'g',  la, A_sysopt3_2 ( : ,  1 )  ,  ' b ' ) 

title ({' \omega=0 . 0004 ,  \mu_l=l/60,  \mu_2=l . 5 \mu_l ,  \mu_3=l . 5\mu_2, . . . 

\nu_0=0 . 0005,  \nu_l=0.001  to  0.01,  \nu_2=l . 5\nu_l ,  \nu_3=l . 5\nu_2 ' ; . . . 
'AVAILABILITY  VS.  ARRIVAL  RATE ' ; ' \nu_l  =  0.001'}) 
legend (' optimal  1  (avail  1)',  'optimal  2  (avail  1)', 'optimal  3  (avail  1)',... 

'optimal  1  (avail  2) ',  'optimal  2  (avail  2) ',  'optimal  3  (avail  2) ') 
xlabel (' Arrival  Rate  (\lambda  1/sec.)'),  ylabel (' Availiability ' ) 
xlim ( [1/300,  1/50] ) 
subplot (1,2,2) 

plot (la, A_sysoptl  (  :  , 2)  ,  ' - . r ' ,  la, A_sysopt2  (  : , 2 ) ,  ' -xg ' ,  la,  A_sysopt 3  (  : , 2 ) ,  .  .  . 

'-ob',  la, A_sysopt 1_2  ( : , 2 ) ,  ' r ' , la, A_sysopt2_2  ( : , 2 ) ,  'g',  la,... 

A_sysopt3_2 ( : , 2 ) , ' b ' ) 

8.4  Expected  response  time 
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state 


=  [  0  0 

0  4 

0  4 

1  1 

1  1 

1  1 

1  2 

1  2 

1  3 

2  1 

2  1 

2  1 

2  2 

2  2 

2  3 

3  1 

3  1 

3  1 

3  2 

3  2 

3  3 

1  2 

1  2 

2  2 

3  2 

1  3 

1  3 

2  3 

2  3 

2  2 

3  3 

3  3 

3  2 

0  4 

1  4 

2  4 

3  4 


0 

0 

4 

0 

4 
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1 
1 
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2 
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5 
5 
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3 
5 
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5 
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%%  parameters 


0 

0 

0 

0 

0 

4 

0 

4 

3 
0 

4 
4 
1 
4 

3 
1 

4 
4 
1 
4 

3 
0 

4 
1 
1 

5 
5 
5 
5 

4 

5 
5 
4 
4 
4 
4 

4]  ; 


la=linspace  (1/300,  0.02,  30); 
nul  =  linspace  (1/1000,  1/100,  2); 
nu0=0.0005;  nu2=1.5*nul;  nu3=1.5*nu2; 
mul=l/60;  mu2=1.5*mul;  mu3=1.5*mu2; 
w=0 . 0004 ; 


for  i=l : length (la) 

for  j=l : length (nul ) 

setpara(mul,  mu2,  mu3,  nuO,  nul(j),  nu2(j),  nu3(j),  la(i),  w)  ; 
failure  =  [nuO  nuO  nuO  nul  nul  nul  nul  nul  nul  nu2  nu2  nu2  nu2  nu2 . . . 
nu2  nu3  nu3  nu3  nu3  nu3  nu3  nul  nul  nu2  nu3  nul  nul  nu2  nu2  nu2 . . . 
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nu3  nu3  nu3  0  0  0  0] ; 


Q  =Qnewl;  %  greedy  case 

Q (:, length (Q) )  =  [ones ( 1 , length (Q) )] ; 
p_ss  =  [ zeros ( 1 ,  (length (Q) -1 )) ,  l]*inv(Q); 

T=P_ss*state ( : , 1 ) ;  %  expected  number  of  targets 

R (i, j ) =T/ (  (3-T) *la  (i) ) ;  %  expected  response  time 

Q2=Qnew2;  %  conservative  case 

Q2 (:, length (Q2 ) )  =  [ones (1, length (Q2) )] ; 
p_ss2  =  [ zeros ( 1 ,  (length (Q2 ) -1 )) ,  l]*inv(Q2); 

T2=p_ss2*state ( : , 1 ) ;  %  expected  number  of  targets 

R2 (i, j ) =T2/ ( (3-T2 ) *la (i) ) ;  %  expected  response  time 

Qag=Qnewag;  %  aggressive  case 

Qag (: , length (Qag) )  =  [ones  (1, length (Qag) )] ; 
p_ssag  =  [zeros (1,  (length (Qag) -1 ))  ,  l]*inv(Qag); 

Tag=p_ssag*state ( : , 1 ) ;  %  expected  number  of  targets 

Rag (i, j ) =Tag/ ( (3-Tag) *la (i) ) ;  %  expected  response  time 

Qopt l=Qnewopt 1 ;  %  optimal  1 

Qopt 1 (:, length (Qoptl ) )  =  [ones ( 1 , length (Qopt 1 ))] ; 
p_ss_optl  =  [zeros (1, (length (Qoptl) -1) ) ,  1 ] *inv (Qopt 1 ) ; 

Topt l=p_ss_opt Instate (: ,  1 )  ;  %  expected  number  of  targets 

Roptl (i, j) =Toptl/ ( (3-Toptl) *la (i) ) ;  %  expected  response  time 

Qopt2=Qnewopt2 ;  %  optimal  2 

Qopt2 (:, length (Qopt2 ) )  =  [ ones ( 1 , length (Qopt2 ))] ; 
p_ss_opt2  =  [zeros (1,  (length (Qopt2) -1) ) ,  l]*inv(Qopt2); 
Topt2=p_ss_opt2 estate (: , 1 ) ;  %  expected  number  of  targets 

Ropt2  (i, j ) =Topt2/ ( (3-Topt2 ) *la (i) ) ;  %  expected  response  time 

Qopt3=Qnewopt3;  %  optimal  3 

Qopt 3 (:, length (Qopt3 ) )  =  [ ones ( 1 , length (Qopt 3 ))] ; 
p_ss_opt3  =  [ zeros ( 1 ,  (length (Qopt3) -1 )) ,  l]*inv(Qopt3); 

Topt 3=p_ss_opt3*state ( : r  1 ) ;  %  expected  number  of  targets 

Ropt3 (i, j) =Topt3/ ( (3-Topt3) *la (i) ) ;  %  expected  response  time 

end, end 

response  time  under  the  breakdown  condition  of  servers 
p=10;  %  assume  10  customers  have  completed  the  service. 

R= (R*p+l/w*ones (30, 2) ) /p; 

R2= (R2*p+l/w*ones (30, 2) ) /p; 

Rag= (Rag*p+l/w*ones (30,2) ) /p; 

Roptl= (Ropt l*p+l/w*ones (30, 2) ) /p; 

Ropt2= (Ropt2*p+l/w*ones (30, 2) ) /p; 

Ropt3= (Ropt3*p+l/w*ones (30, 2) ) /p; 


%%  response  time  vs.  arrival  rate 
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figure  ( 1 ) 
subplot  (1,2, 1 ) 

plot (la, R2 ( : , 1 ) ,  '-o',  la,R(:,l),  ' -x 1 ,  la,Rag(:,l),  ' — ',  la, Ropt 1 ( : , 1 ) , . . . 

la, Ropt2 ( :  ,  1 )  ,  '-+',  la,  Ropt3(:,l)) 

title ({' \omega=0 . 0004 ,  \mu_l=l/60,  \mu_2=l . 5\mu_l ,  \mu_3=l . 5\mu_2 , . . . 

\nu_0  =  0 . 0005,  \nu_l=0.001  to  0.01,  \nu_2  =  l . 5\nu_l ,  \nu_3=l . 5\nu_2 ' ; . . . 
'RESPONSE  TIME  VS.  FAILURE  RATE ' ; 1 \nu_l  =  0.001’}) 
ylabel (' Response  Time  (sec.)'),  xlim( [1/300,  1/50]) 

subplot  (1,2,2) 

plot (la, R2 ( : , 2 ) ,  '-o',  la,R(:,2),  ' -x ' ,  la,Rag(:,2),  ' — ',  la, Ropt 1 ( : , 2 ) , . . . 

la, Ropt2 ( : , 2)  ,  '-+',la,  Ropt3(:,2)) 

legend ( ' conservative ' , ' greedy=optimal  4 ' , ' aggressive ' , ' optimal  1 ' , . . . 

'optimal  2', 'optimal  3') 

title (' \nu_l  =  0.02'),  xlabel (' Arrival  Rate  (\lambda  1/sec.)') 
ylabel (' Response  Time  (sec.)'),  xlim( [1/300,  1/50]) 

8.5  Linear  Programming 


for  k  =  1 : length (state) 

1 (k)  =  state (k,l);  %queue  length  at  each  state 

end 

%%  parameters 

la=l/100;  %  arrival  rate 

w=0.0004;  %  overhaul  rate 

nul=l/1000;  nu0=0.0005;  nu2=1.5*nul;  nu3=1.5*nu2;  %failure  rate 

mul=l/60;  mu2=1.5*mul;  mu3=l . 5*mu2 ;  %service  rate 

rho  =  3* (mul+mu2+mu3+la+nul+nu2+nu3) +w;  %uniform  rate 

beta  =  3;  %  discount  factor 

alpha  =  rho/ (beta+rho) ;  %  discount  factor  for  equivalent  Discrete  Markov 

%  number  of  failed  servers  in  each  state 

m  =  [0120120100120100120101211121221223333]; 
%states  would  lead  to  system  failure 

n  =  [0  01001000001000001000010001011011000  0]  ; 
p  =  [3*nu0  2*nu0  nuO  3*nul  2*nul  nul  2*nul  2*nul  3*nul  3*nu2. . . 

2*nu2  nu2  2*nu2  2*nu2  3*nu2  3*nu3  2*nu3  nu3  2*nu3  2*nu3... 

3*nu3  nul  nul  nu2  nu3  2*nul  nul  2*nu2  nu2  nu2  2*nu3  nu3  nu3  0  0  0  0]; 

step_cost  =  m;  %cost  of  failed  servers 
%  step_cost  =  n/w  ;  %cost  of  replenish 
%  step_cost  =  p;  %cost  of  failure  rate 

%  step_cost  =  1;  %cost  of  queue  length 

setpara(mul,  mu2,  mu3,  nuO,  nul,  nu2,  nu3,  la,  w) ; 

Q1  =  Qmatrix ( 1 , 0 , 0 ) ; 

Q2  =  Qmatrix ( 0 ,  1 ,  0 )  ; 

Q3  =  Qmatrix ( 0 ,  0 ,  1 )  ; 

Pl=eye (37) +Ql/rho;  P2=eye (37) +Q2/rho;  P3=eye  ( 37 ) +Q3/rho;  %P  matrices 
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%%  linear  programming 

f  =  -ones (37,1);  %length  of  f  must  be  same  as  number  of  columns  of  A. 


A  =  -alpha* [Pl-eye (37 ) /alpha;  P2-eye ( 37 ) /alpha; 
P3 (1, 1) -1/alpha  P3(l,2:37); 

P3 (3, 1:2)  P3 (3, 3) -1/alpha  P3(3,4:37); 

P3  ( 5 , 1 :  4 )  P3 (5, 5) -1/alpha  P3(5,6:37); 

P3  (  6, 1 :  5 )  P3 (6, 6) -1/alpha  P3(6,7:37); 

P3 ( 8 , 1 : 7 )  P3 (8, 8) -1/alpha  P3(8,9:37); 

P3  (  9, 1 :  8 )  P3 (9, 9) -1/alpha  P3(9, 10:37); 

P3 ( 10, 1 : 9)  P3 (10, 10) -1/alpha  P3  (10, 11 : 37) ; 
P3  ( 11 , 1 : 10 )  P3 (11, 11) -1/alpha  P3(ll, 12:37) 
P3  ( 12 , 1 : 11 )  P3 (12, 12) -1/alpha  P3(12, 13:37) 
P3  ( 15, 1 : 14 )  P3 (15, 15) -1/alpha  P3(15, 16:37) 
P3  ( 1 6, 1 : 15 )  P3 (16, 16) -1/alpha  P3(16, 17:37) 
P3  ( 17 , 1 : 1 6 )  P3 (17, 17) -1/alpha  P3(17, 18:37) 
P3  ( 18, 1 : 17 )  P3 (18, 18) -1/alpha  P3(18, 19:37) 
P3  (21, 1:20)  P3 (21, 21) -1/alpha  P3(21, 22:37) 
P3  (23,  1:22)  P3 (23, 23 ) -1/alpha  P3  (23, 24:37) 
P3  (2  6,  1:25)  P3 (26, 26) -1/alpha  P3(26,  27:37) 
P3  (27 , 1:26)  P3 (27, 27) -1/alpha  P3(27, 28:37) 
P3 (29, 1:28)  P3 (2 9 , 2 9 ) -1 /alpha  P3(29, 30:37) 
P3 ( 30 , 1:29)  P3 ( 30 , 30 ) -1 /alpha  P3(30, 31:37) 
P3  ( 32 , 1:31)  P3 (32, 32) -1/alpha  P3(32, 33:37) 
P3 ( 33 ,  1:32)  P3 (33,  33) -1/alpha  P3  (33,  34:37) 
P3  ( 34 ,  1:33)  P3 (34 , 34 ) -1/alpha  P3(34, 35:37) 
P3  ( 35,  1 : 34 )  P3 (35,  35) -1/alpha  P3  (35,  36:37) 
P3 (3  6,  1:35)  P3  (36,  36) -1/alpha  P3(36,  37); 

P3 (37,  1 :  36)  P3  (37 , 37 ) -1 /alpha] ; 


b  =[step_cost  step_cost  step_cost(l)  step_cost(3)  step_cost (5) . . . 

step_cost(6)  step_cost(8)  step_cost(9)  step_cost ( 10 )  step_cost ( 11 ) . . . 
step_cost ( 12 )  step_cost (15)  step_cost (16)  step_cost ( 17 )  step_cost (18 ) . . . 
step_cost (21 )  step_cost (23)  step_cost (26)  step_cost (27 )  step_cost  (2 9)  . .  . 
step_cost (30 )  step_cost (32 )  step_cost (33)  step_cost (34 )  step_cost (35) . . . 
step_cost (36)  step_cost (37 ) ] ' ;  %length  b  must  be  =  number  of  rows  of  A. 

lb  =  zeros (37,1); 

[x, fval, exitf lag, output , lambda]  =  linprog ( f , A, b, [ ] , [ ] , lb) ; 

%%  comparison 
for  i=l : 37 

vl (i) =step_cost (i) talpha* (PI (i, : ) *x) ; 
v2 (i) =step_cost (i) talpha* (P2 (i, : ) *x) ; 
v3  (i) =step_cost (i) talpha* (P3 (i,  : ) *x) ; 
al ( i ) = v 1 (i) -x  (i)  ; 
a2 ( i ) =v2 ( i ) -x  ( i ) ; 
a3 ( i ) =v3 ( i ) -x  ( i ) ; 

end 


%display  optimal  cost  and  cost  from  the  three  control  sets 
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disp  (  1  opt_value  vl  v2  v3 f  ) 

disp ( [x  vl 1  v2 1  v3 * ] ) 

%  determine  the  control  policy  for  each  state  by 
%  checking  out  the  smallest  value  among  al  a2  a3 
disp  (  1  control  policy:  '  ) 
disp ( ' ul  u2  u3 1 ) 
for  i=l : 37 

if  (abs  (al  (i)  '  )  >  abs(a3(i)'))  &  (abs(a2(i)')  >  abs(a3(i)’)) 
disp  ( 1  0  0  1') 

elseif  (abs(al(i)’)  >  abs(a2(i)'))  &  (abs(a3(i)')  >  abs(a2(i)')) 
disp  ( 1  0  1  O') 

elseif  (abs(a2(i)')  >  abs(al(i)'))  &  (abs(a3(i)')  >  abs(al(i)')) 
disp  (  '  1  0  O') 

else  disp  ( 1  0  0  O') 

end 

end 
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Red:  arrival  and  departure 

Green:  failure 
Blue:  departure  only 

Purple:  overhaul 


9  MATLAB  codes  and  simulations  of  the 
airborne  Network 


9.1  Closed  queuing  network  (reducible) 


Figure  9.1:  Top  level  Simulink  block  diagram  of  simulation  model  1 
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9.1.1  Stateflow  diagram  of  the  system  state  (reducible) 


Figure  9.2:  System  Stateflow  diagram  of  simulation  model  1 
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9.1.2  Stopping  rule 


Do  Not  Delete 


Subsystem  Configuration 


Figure  9.3:  Stopping  rule  of  simulation  model  1 


9.1.3  Computing  total  number  of  targets  in  the  system 

Subsystem  Configuration 


Din 

Figure  9.4:  Computing  total  number  of  targets 


9.1.4  Feedback  subsystem 


Figure  9.5:  Lower  portion  subsystem 
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9.1.5  Identify  system’s  next  state 


|  Do  Not  Delete  | 

Subsystem  Configuration 

Identify  next  state 


Figure  9.6:  Identify  next  state 


The  following  is  the  embedded  MATLAB  code  for  the  next  state  subsystem: 


function  next  =  sysl (state, dl, d2, d3) 
persistent  last_state; 

%  states :  (dl ,  d2 ,  d3 )  ,  0  represents  sensor  up,  1  represents  sensor  down. 
%  0  —  000 

%  1  —  001 

%  2  —  010 

%  3  —  Oil 

%  4  —  100 

%  5  —  101 

%  6  —  110 

%  7  —  111 

%  initialize  output 
next=0 ; 
switch  state 
case  0 

xs= [dl; d2; d3] ; 
states=  [4  2  1 ] ; 

if  xs ( 1 ) — — 1  &&  xs (2) ==1  &&  xs (3 ) ==1 
next=7 ; 

elseif  xs ( 1 ) ==1  &&  xs(2)==l  &&  xs(3)==0 
next=6; 

elseif  xs  ( 1 ) ==1  &&  xs(2)==0  &&  xs(3)==l 
next=5 ; 
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elseif  xs ( 1 ) ==0  &&  xs ( 2 ) ==1  &&  xs(3)==l 
next=3 ; 

elseif  xs  ( 1 ) ==0  &&  xs(2)==0  &&  xs(3)==0 
next=0 ; 

else  [x  ind] =max (xs ) ; 

next=states (ind) ; 
end 

case  1 

xs=[d2;  dl;  -inf]; 
states^ [3  5  0  ]  ; 
if  xs  ( 1 )  ==0  &  Sc  xs  ( 2 )  ==0 
next=l ; 

elseif  xs ( 1 ) ==1  &&  xs ( 2 ) ==1 
next=7 ; 

else  [x  ind] =max (xs ) ; 
next=states (ind) ; 

end 

case  2 

xs=[d3;  dl;  -inf]; 
states^ [3  6  0 ] ; 
if  xs(l)==0  Sc  Sc  xs  (2)  ==0 
next=2 ; 

elseif  xs ( 1 ) ==1  &&  xs(2)==l 
next=7 ; 

else  [x  ind] =max (xs ) ; 
next=states (ind)  ; 
end 

case  4 

xs=[d3;  d2;  -inf]; 
states= [5  6  0 ] ; 
if  xs  ( 1 )  ==0  Sc  Sc  xs  ( 2 )  ==0 
next=4 ; 

elseif  xs ( 1 ) ==1  &&  xs ( 2 ) ==1 
next=7 ; 

else  [x  ind] =max (xs )  ; 

next=states (ind) ; 
end 

case  3 

if  dl==l 

next=7;  end 

case  5 

if  d2==l 

next=7;  end 

case  6 

if  d3==l 

next=7;  end 

case  7 
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end 


next=7 ; 


9.1.6  Service  rate  generator 


service  rate  generator 


Subsystem  Configuration 


Figure  9.7:  Service  rate  generator 


function  [ rl ,  r2 ,  r3 ,  cl ,  c2 ,  c3 ]  =  system (Q, mul , mu2 , mu3 , FI , F2 , F3 ) 

%  Q  =  queue  length 

%  FI,  F2 ,  F3 :  state  of  sensors.  0 — >UP,  1 — >DOWN 
%  cl,c2,c3:  enable  gate  signal 

%  This  script  contains  four  control  policies.  Run  the  simulation 
%  at  each  control  policy  and  plot  the  results  together. 


%%GREEDY 

%  if  F1==0  &&  F2==0  &&  F3==0 

%  rl  =  mu3;  r2=mu3;  r3=mu3; 

%  cl=l;  c2=0;  c3=0; 

%  elseif  F1==0  &&  F2==0  &&  F3==l 

%  rl  =  mu2;  r2=mu2;  r3=0; 

%  cl=l;  c2=0;  c3=0; 

%  elseif  F1==0  &&  F2==l  &&  F3==0 

%  rl  =  mu2;  r2=0;  r3=mu2; 

%  cl=l;  c3=0;  c2=0; 

%  elseif  Fl==l  &&  F2==0  &&  F3==0 

%  r2  =mu2;  rl=0;  r3=mu2; 

%  c2=l;  c3=0;  cl=0; 

%  elseif  F1==0  &&  F2==l  &&  F3==l 

%  rl  =  mul;  r2=0;  r3=0; 

%  cl=l;  c2=0;  c3=0; 

%  elseif  Fl==l  &&  F2==0  &&  F3==l 

%  r2  =  mul;  rl=0;  r3=0; 

%  c2=l;  cl=0;  c3=0; 

%  elseif  Fl==l  &&  F2==l  &&  F3==0 

%  r3  =  mul;  rl=0;  r2=0; 
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%  c3=l;  cl=0;  c2=0; 

%  else  rl=0;  r2=0;  r3=0; 

%  cl=0;  c2  =  0;  c3=0; 

%  end 

0.0. _ 

o  o 


%  %%  CONSERVATIVE 
%  if  F1==0  &&  F2==0  &&  F3==0 
%  rl  =  mul;  r2  =  mul;  r3  =mul; 

%  cl=l;  c2=l;  c3=l; 

%  elseif  F1==0  &&  F2==0  &&  F3==l 

%  rl  =  mul;  r2=  mul;  r3  =  mul; 

%  cl=l;  c2=l;  c3=0; 

%  elseif  F1==0  &&  F2==l  &&  F3==0 

%  rl  =  mul;  r3  =  mul;  r2  =mul; 

%  cl=l;  c3=l;  c2=0; 

%  elseif  Fl==l  &&  F2==0  &&  F3==0 

%  r2  =mul;  r3  =  mul;  rl  =  mul; 

%  c2=l;  c3=l;  cl=0; 

%  elseif  F1==0  &&  F2==l  &&  F3==l 

%  rl  =  mul;  r2  =  mul;  r3  =  mul; 

%  cl=l;  c2=0;  c3=0; 

%  elseif  Fl==l  &&  F2==0  &&  F3==l 

%  r2  =  mul;  rl  =  mul;  r3  =  mul; 

%  c2=l;  cl=0;  c3=0; 

%  elseif  Fl==l  &&  F2==l  &&  F3==0 

%  r3  =  mul;  rl  =  mul;  r2  =  mul; 

%  c3=l;  cl=0;  c2=0; 

%  else  rl=0;  r2=0;  r3=0; 

%  cl=0;  c2=0;  c3=0; 

%  end 

0,0, _ 

o  o 


%  AGGRESSIVE 

%  if  F1==0  &&  F2==0  &&  F3==0  &&  Q==0 


Q. 

O 

rl 

=  mu  3; 

r2=mu3;  r3=mu3; 

Q. 

O 

c  1  =  1  ; 

c2=l;  c3=l; 

Q. 

O 

elseif 

F1==0 

&&  F2==0  &&  F3==0 

&& 

Q==l 

Q. 

O 

rl 

=  mu  3; 

r2=mu3;  r3=mu3; 

Q. 

O 

c  1  =  1  ; 

o 

N) 

II 

O 

O 

CO 

II 

o 

O. 

O 

elseif 

F1==0 

&&  F2==0  &&  F3==0 

&& 

Q==2 

Q. 

O 

rl 

=  mu2 ; 

r2=mu2;  r3=mul; 

Q. 

O 

c  1  =  1  ; 

c2=0;  c3=l; 

O. 

O 

elseif 

F1==0 

&&  F2==0  &&  F3==0 

&& 

Q>3 

Q. 

O 

rl 

=  mul; 

r2=mul;  r3=mul; 

O. 

O 

Q. 

c  1  =  1  ; 

c2=l;  c3=l; 

O 

Q. 

O 

elseif 

F1==0 

&&  F2==0  &&  F3==l 

&& 

Q==0 

Q. 

O 

rl 

=  mu2 ; 

r2=mu2;  r3=0; 

O. 

O 

c  1  =  1  ; 

o 

K) 

II 

I—4 

O 

CO 

II 

o 

Q. 

O 

elseif 

F1==0 

&&  F2==0  &&  F3==l 

&& 

Q==l 
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rl 

=  mu2 ; 

r2=mu2 

o 

II 

cn 

Ci 

cl  =  l ; 

c2  =  0 ; 

c3=0 ; 

elseif 

M ‘ 

II 

II 

O 

II 

II 

c\j 

Cm 

o 

CO 

II 

II 

h-4 

&& 

Q==2 

rl 

=  mul ; 

r2=mul 

;  r3=mul; 

c  1  =  1 ; 

c2  =  l ; 

c3=0 ; 

elseif 

I—4 

II 

II 

o 

II 

II 

CN 

Cm 

<-3 

<-3 

0  &&  F3==l 

&& 

Q>3 

rl 

=  mul; 

r2=mul 

;  r3=mul; 

c  1  =  1 ; 

c2  =  l ; 

o 

II 

CO 

o 

elseif 

I—4 

II 

II 

o 

II 

II 

CM 

Pm 

<-3 

1  &&  F3==0 

&& 

Q==0 

rl 

=  mu2 ; 

o 

II 

CM 

Cl 

r3=mu2 ; 

c  1  =  1 ; 

O 

CO 

II 

o 

c2  =  l  ; 

elseif 

o 

II 

II 

\ — 1 
Cm 

&&  F2== 

1  &&  F3==0 

&& 

Q==l 

rl 

=  mu2 ; 

o 

II 

CM 

Cl 

r3=mu2 ; 

c  1  =  1 ; 

O 

CO 

II 

o 

c2=0 ; 

elseif 

I—4 

II 

II 

o 

II 

II 

c\] 

Cm 

I-4 

CO 

ll 

ll 

o 

&& 

Q==2 

rl 

=  mul; 

ro 

II 

o 

r3=mul ; 

c  1  =  1 ; 

c3=l ; 

c2=0 ; 

elseif 

I—4 

II 

II 

o 

II 

II 

c\j 

Pm 

Ob 

1  &&  F3==0 

&& 

Q>3 

rl 

=  mul; 

r2=0 ; 

r3=mul ; 

\ — 1 

II 

\ — 1 
o 

c3=l ; 

c2=0 ; 

elseif 

I—4 

II 

II 

h-4 

&&  F2== 

o 

CO 

II 

II 

o 

&& 

Q==0 

r2  =mu2;  rl=0;  r3=mu2; 
c2=l;  c3=l;  cl=0; 

elseif  Fl==l  &&  F2==0  &&  F3==0  &&  Q==l 
r2  =mu2;  r 1=0 ;  r3=mu2; 
c2=l;  c3=0;  cl=0; 

elseif  Fl==l  &&  F2==0  &&  F3==0  &&  Q==2 
r2  =mul;  rl=0;  r3=mul; 
c2=l;  c3=l;  cl=0; 

elseif  Fl==l  &&  F2==0  &&  F3==0  &&  Q>3 
r2  =mul;  rl=0;  r3=mul; 
c2=l;  c3=l;  cl=0; 

elseif  F1==0  &&  F2==l  &&  F3==l 
rl  =  mul;  r2=0;  r3=0; 
cl=l;  c2=0;  c3=0; 

elseif  Fl==l  &&  F2==0  &&  F3==l 
r2  =  mul;  r 1 = 0 ;  r3=0; 
c2=l;  cl=0;  c3=0; 

elseif  Fl==l  &&  F2==l  &&  F3==0 
r3  =  mul;  rl=0;  r2=0; 
c3=l;  cl=0;  c2=0; 

else  rl=0;  r2=0;  r3=0; 

cl=0;  c2=0;  c3=0; 

end 


Optimal  1 
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if  F1==0  ScSc  F2==0  &&  F3==0  &&  Q==0 
rl  =  mu2;  r2=mu2;  r3=mu2; 
cl=l;  c2=l;  c3=l; 

elseif  F1==0  &&  F2==0  &&  F3==0  &&  Q==l 
rl  =  mu2;  r2=mu2;  r3=mul; 
cl=l;  c2=0 ;  c3=l; 

elseif  F1==0  &&  F2==0  &&  F3==0  &&  Q==2 
rl  =  mu2;  r2=mu2;  r3=mul; 
cl=l;  c2  =  0;  c3=l; 

elseif  F1==0  &&  F2==0  &&  F3==0  &&  Q>3 
rl  =  mu2;  r2=mu2;  r3=mul; 
cl=l;  c2  =  0;  c3=l; 

elseif  F1==0  &&  F2==0  &&  F3==l  &&  Q==0 
rl  =  mu2;  r2  =mu2;  r3  =  mul; 
cl=l;  c2  =  0;  c3=0; 

elseif  F1==0  &&  F2==0  &&  F3==l  &&  Q==l 
rl  =  mu2;  r2  =mu2;  r3  =  mul; 
cl=l;  c2=0;  c3=0; 

elseif  F1==0  &&  F2==0  &&  F3==l  &&  Q>2 
rl  =  mul;  r2  =mul;  r3  =  mul; 
cl=l;  c2=l;  c3=0; 

elseif  F1==0  &&  F2==l  &&  F3==0  &&  Q==0 
rl  =  mu2;  r2=0;  r3=mu2; 
cl=l;  c3  =  0;  c2  =  0; 

elseif  F1==0  &&  F2==l  &&  F3==0  &&  Q==l 
rl  =  mu2;  r2  =  0;  r3=mu2; 
cl=l;  c3  =  0;  c2  =  0; 

elseif  F1==0  &&  F2==l  &&  F3==0  &&  Q>2 
rl  =  mul;  r2=0;  r3=mul; 
cl=l;  c3=l;  c2=0; 

elseif  Fl==l  &&  F2==0  &&  F3==0  &&  Q==0 
r2  =mu2;  rl  =  0;  r3=mu2; 
c2=l;  c3  =  0;  cl=0; 

elseif  Fl==l  &&  F2==0  &&  F3==0  &&  Q==l 
r2  =mu2;  rl  =  0;  r3=mu2; 
c2=l;  c3  =  0;  cl=0; 

elseif  Fl==l  &&  F2==0  &&  F3==0  &&  Q>2 
r2  =mu 1 ;  rl  =  0;  r3=mul; 
c2=l;  c3=l;  cl  =  0; 

elseif  F1==0  &&  F2==l  &&  F3==l 
rl  =  mul;  r2=0;  r3  =  0; 
cl=l;  c2=0;  c3=0; 

elseif  Fl==l  &&  F2==0  &&  F3==l 
r2  =  mul;  rl=0;  r3  =  0; 
c2=l;  cl  =  0;  c3=0; 

elseif  Fl==l  &&  F2==l  &&  F3==0 
r3  =  mul;  rl=0;  r2  =  0; 
c3=l;  cl=0;  c2  =  0; 
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else  rl=0;  r2=0;  r3=0; 
cl=0;  c2=0;  c3=0; 

end 


9.1.7  New  allocation  for  return  targets 


Do  Not  Delete 
Subsystem  Configuration 


MATLAB  Functionl 

Figure  9.8:  New  allocation  for  return  targets 


function  [ccl, cc2, cc3]  =  sys2 ( state, next , pi , p2 , p3 , restart ) 

%  pl,p2,p3  are  the  ejected  entities  from  servers 
%  Q:  upper  system  queue  length,  include  the  ones  in  service 
%  state:  current  state  of  the  system 
%  next:  next  state  of  the  system 

%  ccl,cc2,cc3:  signals  determine  which  server  serves  a  restarted  entity 

o,  _ 
o 

if  restart  >1 

if  state  ==0  &&  next==l  &&  p3==l 
ccl=l;  cc2=0;  cc3=0; 
elseif  state==0  &&  next==2  &&  p2==l 
ccl=l;  cc2=0;  cc3=0; 
elseif  state==0  &&  next==4  &&  pl==l 
ccl=0;  cc2=l;  cc3=0; 

elseif  state  ==1  &&  next==3  &&  p2==l 
ccl=l;  cc2=0;  cc3=0; 
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elseif  state  ==1  &&  next==5  &&  pl==l 
ccl=0;  cc2=l;  cc3=0; 

elseif  st ate==2  &&  next==3  &&  p3==l 
ccl=l;  cc2=0;  cc3=0; 

elseif  st ate==2  &  &  next==6  &  &  pl==l 
ccl=0;  cc2=0;  cc3=l; 

elseif  state==4  &  &  next==5  &  &  p3==l 
ccl=0;  cc2=l;  cc3=0; 

elseif  st ate==4  &&  next==6  &&  p2==l 
ccl=0;  cc2=0;  cc3=l; 

elseif  state==3  &&  next==7 

ccl=0;  cc2=0;  cc3=0; 

elseif  state==5  &&  next==7 

ccl=0;  cc2=0;  cc3=0; 

elseif  state==6  &&  next==7 

ccl=0;  cc2=0;  cc3=0; 

else  ccl=0;  cc2=0;  cc3=0; 

end 

else  ccl=0;  cc2=0;  cc3=0; 
end 
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9.1.8 


Sensor-pair  subsystem 


Figure  9.9:  Sensor-pair  subsystem  for  simulation  model  1 
The  system  has  three  identity  sensor-pairs  subsystem. 
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9.1.9  Failure  generation  subsystem 


l 


Figure  9.10:  Failure  generation  subsystem 


All  sensor-pairs  have  same  random  failure  time;  therefore  all  sensor-pairs  subsystems  contain 
the  following  subsystem. 
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9.1.9.1 


Discrete  event  subsystem! 


Do  Not  Delete 

Subsystem  Configuration 

Embedded 
MATLAB  Function 

-KD 

Dout 


Figure  9.11:  Discrete  event  subsysteml 

This  block  is  within  the  failure  generation  subsystem,  follow  code  is  embedded  in  the  block. 
The  code  generates  exponential  failure  time  of  the  sensor-pairs. 

function  muhat  =  fcn(t) 
persistent  oldtime; 

%initialize  output 
muhat=2000 ; 

%initialize  persistent  variable 
if  isempty (oldtime) 
oldtime=0 ; 

end 

%update 
if  (t>oldtime) 
oldtime=t ; 

muhat =  expfit (exprnd (2000,  10,  1)  )  ; 
end 
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9.1.9.2  Discrete  event  subsystem  2 


Embedded 
MATLAB  Function 


Dout 


Figure  9.12:  Discrete  event  subsystem  2 

The  above  block  is  also  within  failure  generation  subsystem,  following  code  is  embedded  in 
the  block  which  assigns  failure  time  to  failure  entity. 

function  t  =  fcn(ul,u2) 
persistent  oldtime; 

%initialize  output 
t=0 ; 

%initialize  persistent  variable 
if  isempty (oldtime) 
oldtime=0 ; 

end 

%update 

if  (ul>oldtime) 
oldtime=ul ; 
end 

if  ul>u2 
t  =  0  ; 

else  t=u2-oldtime; 
end 
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9.1.9.3  Discrete  event  subsystem  3 


Din 

Figure  9.13:  Discrete  event  subsystem  3 


If  there  is  an  entity  exits  the  server  in  the  failure  generation  subsystem  of  a  sensor-pair,  then 
that  sensor-pair  fails.  Otherwise,  the  sensor-pair  is  still  up. 
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9.2  Irreducible  case  with  system  replenishment  and  with  or 
without  channel  fading 
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9.2.1  Stateflow  diagram  of  the  system  (irreducible) 

The  stateflow  diagram  is  similar  to  fig  C.2,  except  it  has  a  transition  from  final  failure  state  to  initial 
state. 


Figure  9.15:  System  state-flow  diagram  for  simulation  model  2  and  3. 
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9.2.2  System  up/down  stateflow  diagram 


Figure  9.16:  State-flow  diagram  of  network  status. 


9.2.3  Stopping  rule  (irreducible) 


Do  Not  Delete 


Subsystem  Configuration 


► 


Compare  Data  Type  Conversion  Stop  Simulation 
To  Constant 


Figure  9.17:  Stopping  rule  for  simulation  model  2  and  3. 


co — ►  ==6 


9.2.4  Repairing  rule 


|  Do  Not  Delete  j 
Subsystem  Configuration 


MATLAB  Function 


Figure  9.18:  Repairing  rule  of  simulation  model  2  and  3. 
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9.2.5  Sensor-pair  subsystem  (irreducible)  without  channel  fading 
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9.2.6  Failure  generation  subsystem  (irreducible) 


Figure  9.20:  Failure  generation  subsystem  for  simulation  model  2  and  3. 
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9.2.6. 1  Discrete  event  subsystem  2 
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Dout 


Figure  9.21:  This  block  is  within  Fig.  C.16,  the  function  of  this  block  generates  failure  time. 


function  t  =  f cn (ul , u2 , elapsed, past ime) 
persistent  oldtime; 

%initialize  output 
t=0 ; 

%initialize  persistent  variable 
if  isempty (oldtime) 
oldtime=0 ; 

end 

%update 

if  (ul>oldtime) 
oldtime=ul ; 
end 

T  =  elapsed+pastime; 
if  u2  <  ul-T 
t  =  0  ; 

else  t=u2-ul+T; 
end 
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9.2.6.2  Discrete  event  subsystem  3 
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Figure  9.22:  Discrete  event  subsystem  3 


1  function  dout=fcn(d) 

2 

3  if  d==0 

4  dout=l; 

5  else  dout=0; 

6  end 
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9.3  Irreducible  case  with  system  replenishment  and  intermit¬ 
tent  channel  fading 

This  model  is  similar  to  model  2  except  intermittent  channel  fading  exists;  the  top  level  block 
diagram  is  the  same  as  model  2’s.  This  section  illustrates  the  additional  blocks  to  the  sensor-pair 
subsystem  shown  in  the  following  figure. 


Figure  9.23:  Sensor-pair  subsystem  of  simulation  model  3 


Figure  C.  18  shows  how  to  setup  the  temporary  outage  and  how  to  recovery  the  loss  data  for  the 
sensor  pairs  and  figure  C.19  determines  the  recovery  rate. 
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Figure  9.24:  Outage  and  recovery  subsystem 
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Figure  9.25:  Discrete-event  subsystem  for  recovery  rate 


function  t  =  fen (n, gamma) 
if  n==0 
t  =  0  ; 

else  t=gamma; 
end 


9.4  Code  for  running  the  simulations  repeatedly 

9.4.1  Initial  seed  value  generator: 


%  Generate  a  vector  of  large  odd  numbers, 
newseed  =  (30000  :  2  :  79999); 

%  Randomly  permute  the  numbers  to  avoid  always  using  the  same  set  of  seeds, 
perm  =  randperm ( length (newseed) ) ; 

paramname  =  { ' initialseed '  ,  'seed'};  %  Parameter  names  to  consider 
np  =  length (paramname) ; 
idx  =  1; 

for  j  j  =  l : np 

%  Find  blocks  that  have  the  parameter  with  a  numerical  value 
blocks  =  f ind_system (bdroot , 'RegExp ' , ' on ' , ' LookUnderMasks ' , ' all ' , . . . 
paramname {  j  j  } ,  ' \d+ ' ) ; 
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%  Replace  initial  seed  parameter  values  with  numbers  from  the  vector, 
for  kk  =  1 : length (blocks ) 

newseedparamvalue  =  num2str (newseed (perm (idx) ) ) ; 
idx  =  idx  +  1; 

set_param (blocks { kk } , paramname { j  j } ,  newseedparamvalue) ; 
disp ([ 1  Setting  parameter  to  '  newseedparamvalue  '  in  '  ... 

strrep (blocks { kk }, char ( 10 ),  ’  ') 

end 

end 

9.4.2  Run  simulation  repeated  to  gather  results  and  varying  a  parameter 

9.4.2.1  1.  Mean  time  to  network  expiration 


%  execute  the  following  two  comments  in  MATLAB,  resave  the  model,  then 
%  start  the  simulation. 

%  set_param ( ' closetestl ' ,  ' PostLoadFcn ' ,  '  ' ) 

%  set_param ( ' closetestl ' , ' CloseFcn ' , ' ' ) 

clear  global  all; 

load_system (' closetestl 1 ) ;  %  Load  system 
nruns  =30;  %  Number  of  simulation  runs 

%initialization 

w  =  zeros (nruns, 1) ;  %MTTNE 

d  =  zeros (nruns, 1) ;  %#  of  entity  departed 

opts  =  simset ( ' SrcWorkspace ' ,  ' Current '  ,  ' DstWorkspace '  ,  ' Current ' )  ; 
h  =  waitbar ( 0 Running  simulations.  Please  wait... Be  patient'); 
seedarrayl  =  ceil (rand (nruns , 1 )* 99999) *2+1 ; 
seedarray2  =  ceil (rand (nruns , 1 )* 99999) *2+1 ; 

%  Vary  the  mean  arrival  rate. 

la  =  linspace (1/300, 0 . 02,  4 )  ;  %  Values  of  mean  arrival  rate 
wavg  =  zeros (length (la) , 1) ; 
davg  =  zeros (length (la) , 1) ; 
for  i  =  1: length (la) 

%  m  is  a  parameter  in  the  random  number  block,  so  changing  m  changes  the 
%  mean  arrival  rate  in  the  simulation, 
m  =  la ( i )  ; 

disp ([' Simulating  with  mean  arrival  rate='  num2str (m) ] ) ; 

%  Replicate  for  each  value  of  m 
for  k  =  1 : nruns 

waitbar ( k / nruns , h) ; 
seedvaluel  =  seedarrayl (k) ; 
seedvalue2  =  seedarray2  (k) ; 

sim (' closetestl ',[], opts ) ;  %  Run  simulation, 
w (k)  =  MTTNE ; 
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d(k)  =  entity; 


end 

wavg(i)  =  mean (w) ;  %  Average  for  fixed  value  of  m 

davg(i)  =  mean (d) ; 

end 

close (h) ; 

disp ( 1  average  #  of  entity  departed:  ') 
disp (davg) 
disp (wavg) 

%%  plot  avg.  MTTNE 
figure; 

plot (la, wavg,  '-o'); 

title ('Mean  Time  to  Network  Failure  vs.  Arrival  Rate') 
xlabel (' Arrival  Rate'),  ylabel (' MTTNE  (sec)') 

9.4.2.2  Response  time  and  availability 

The  following  code  is  for  design  2,  may  change  the  model  name  for  design  3. 

%  execute  the  following  two  comments  in  MATLAB,  resave  the  model,  then 
%  start  the  simulation. 

%set_param ( ' closed_replenishment3 ' ,  ' PostLoadFcn ' ,  '  ' ) 

%set_param ( ' closed_replenishment3 ' , ' CloseFcn ' , ' ' ) 

clear  global  all; 

load_system ( ' closed_replenishment3 ' ) ;  %  Load  system 
nruns  =30;  %  Number  of  simulation  run 

%initialization 

R  =  zeros (nruns, 1) ;  %response  time 
A  =  zeros (nruns, 1) ;  %availability 

opts  =  simset ( ' SrcWorkspace ' , ' Current ' , ' DstWorkspace ' , ' Current ' ) ; 
h  =  waitbar ( 0 Running  simulations.  Please  wait...'); 
seedarrayl  =  ceil ( rand (nruns , 1 )* 99999) *2+1 ; 
seedarray2  =  ceil ( rand (nruns , 1 )* 99999) *2+1 ; 
seedarray3  =  ceil ( rand (nruns , 1 )* 99999) *2+1 ; 

la  =  linspace (1/300, 1/50, 4 ) ;  %  Values  of  mean  arrival  rate 
Ravg  =  zeros (length (la) , 1) ; 

Aavg  =  zeros (length (la) , 1) ; 

%  Vary  the  mean  arrival  rate, 
for  i  =  1: length (la) 

%  m  is  a  parameter  in  the  random  number  block,  so  changing  m 
%  changes  the  mean  arrival  rate  in  the  simulation, 
m  =  la ( i )  ; 

disp ([' Simulating  with  mean  arrival  rate='  num2str (m) ] ) ; 

%  Replicate  for  each  value  of  m 
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for  k  =  1  inruns 

wait bar (k/ nruns ,h) ; 
seedvaluel  =  seedarrayl (k) ; 
seedvalue2  =  seedarray2 (k) ; 
seedvalue3  =  seedarray3 (k) ; 

sim ( ' closed_replenishment 3 ' , [ ] ,  opt s ) ;  %  Run  simulation. 

R(k)  =  response_t ime; 

A(k)  =  avail; 

end 

Ravg(i)  =  mean (R) ;  %  Average  for  fixed  value  of  R  and  A 

Aavg(i)  =  mean (A) ; 

end 

close (h) 

%  plot  response  time  and  availability 

figure (1),  plot (la, Ravg, ' -o ' ) ; 

t it le (' Response  Time  vs.  Arrival  Rate') 

xlabel (' Arrival  Rate'),  ylabel (' Response  Time  (sec)') 

figure (2),  plot (la, Aavg, ' -x 1 ) ; 

tit le (' Availability  vs.  Arrival  Rate') 

xlabel (' Arrival  Rate'),  ylabel (' Availability ' ) 
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10  Symbols,  Abbreviations  and  Acronyms 


FDOA  Frequency  Difference  of  Arrival 
MTTNE  Mean  Time  to  Network  Expiration 
TDOA  Time  Difference  of  Arrival 
UAV  Unmanned  Aerial  Vehicle 
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